AN  INVESTIGATION  OF  DIGITAL  SPECTRAL  ANALYSIS 
PROGRAMS  AND  COMPUTER  METHODS  UTILIZED  AT  THE 
NAVAL  POSTGRADUATE  SCHOOL  IN  THE  ANALYSIS  OF 
HIGH  FREQUENCY  RANDOM  SIGNALS 


John  DeMille  McKendrick 


NAVAL  POSTGRADUATE  SCHOOL 

Monterey,  California 


THESIS 


AN  INVESTIGATION  OF  DIGITAL  SPECTRAL  ANALYSIS 
PROGRAMS  AND  COMPUTER  METHODS  UTILIZED  AT  THE 
NAVAL  POSTGRADUATE  SCHOOL  IN  THE  ANALYSIS  OF 
HIGH  FREQUENCY  RANDOM  SIGNALS 

by 

John  DeMille  McKendrick 

Thesis  Advisor:             N 

.  E. 

J.  Boston 

March      1972 


Approved  ^oh.  pubtic  kqJLqjo&z;   dibtxibuXion  untunltzd. 


An  Investigation  of  Digital  Spectral  Analysis  Programs 
and  Computer  Methods  Utilized  at  the  Naval  Postgraduate 
School  in  the  Analysis  of  High  Frequency  Random  Signals 


by 


John  DeMille  McKendrick 
Lieutenant,  United  States  Navy 
B.S.,  United  States  Naval  Academy,  1966 


Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 


MASTER  OF- SCIENCE  IN  OCEANOGRAPHY 


from  the 


NAVAL  POSTGRADUATE  SCHOOL 
March  1972 


ABSTRACT 

The  digitizing  procedure  used  at  the  Naval  Postgraduate 
School  was  investigated  for  possible  sources  of  noise  and  other 
errors.   Signals  of  known  form  were  digitized  through  the 
Analog-to  Digital  Hybrid  computer  system  (Ci  5000/XDS9300)  . 
Similar  signals  were  generated  by  digital  programs  on  the 
IBM  36O/67.   The  resultant  signals  were  analyzed  by  the  com- 
puter programs  UBCFTOR,  which  computed  the  Fourier  coefficients 
of  each  block  of  data,  and  by  UBCSCOR,  which  computed  the 
power  spectra  of  the  signals.   The  power-spectral  plots  of 
the  computer-generated  signals  were  compared  with  the  power- 
spectral  plots  of  digitized  signals.   The  analog-to-digitital 
process  appeared  to  be  relatively  noise  free. 

To  further  test  the  system,  atmospheric  temperature  and 
wind  velocity  signals  were  digitized  and  analyzed  under 
UBCFTOR  and  UBCSCOR.   Plots  of  the  time-varying  spectra  of 
these  signals  compared  favorable  with  results  obtained  at 
other  digitizing  facilities. 
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I.   INTRODUCTION 

A.  PROBLEM 

Random  processes  are  often  recorded  in  a  fluctuating  yol- 
tage  in  analog  form.   However,  it  is  frequently  more  convenient 
to  analyze  the  signal  digitally  (on  large  digital  computers). 
Thus,  investigators  are  faced  with  the  problem  of  converting 
a  continuous  signal  into  discrete  data  samples,  and  with  sub- 
sequent analysis  of  these  digitized  samples.   There  are 
several  steps  in  the  analog-to-digital  conversion  procedure, 
and  in  the  digital  analysis  procedure,  which  allow  for  errors 
and  for  possible  contamination  of  a  signal  with  noise  from 
external,  sources. 

B.  OBJECTIVE 

The  main  goal  of  this  study  was  to  investigate  digitization 
and  analysis  procedures  for  possible  sources  of  noise  which 
may  be  introduced  to  the  true  signals.   The  overall  procedure 
used  at  the  Naval  Postgraduate  School,  from  digitaization  of 
the  actual  signal  to  the  computation  of  the  Power-Spectral- 
Density  (PSD),  is  quite  complex  and  requires  four  computer 
programs.   The  next  objective  was  to  improve  the  routine 
procedures  required  in  this  particular  time  series  analysis 
technique.   The  final  objective  was  to  digitize  actual 
geophysical  signals  and  compare  the  results  with  similar 
analyses  of  these  signals  undertaken  at  the  University  of 
British  Columbia,  by  Boston  in  1970  [Ref.  1]. 
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II.   THEORY 

A.   DIGITAL  REPRESENTATION  OF  CONTINUOUS,   TIME-VARYING  SIGNAL 

Random  geophysical  processes  are  often  studied  by  re- 
cording a  continually  changing  event  as  a  continuous,  fluctuat- 
ing voltage,  which  corresponds  linearly  with  the  original 
process.   The  actual  geophysical  variables  considered  later 
in  this  study  were  small-scale  fluctuations  of  air  temperature, 
wind  velocity,  and  time  derivatives  of  both. 

1.  Analog-to-Digital  Conversion 

A  randomly  fluctuating  voltage  signal  might  look  like 
the  one  in  Figure  1(a).   In  order  to  analyze  the  signal  by 
digital  techniques,  discrete  samples  of  the  fluctuating  voltage 
must  be1  taken,  and  these  are  referred  to  as  sequential  digital 
samples  (vi).   The  requirement  for  sampling  at  equal  intervals 
of  time  is  set  by  the  assumption,  in  most  analyses,  that  they 
are  equal  time  interval  samples.   The  digitized  samples  would 
look  like  Figure  1(b). 

2 .  ,  Digital  Sampling  Theory 

Implied,  in  the  digitization  problem,  is  ttie  question 
of  how  well  the  sequence  represents  the  original  voltage.  In 
turning  to  Sampling  Theory  for  the  answer,  three  hypotheses 
are  made: 

a.  V(t)  is  a  random  variable  defined  for  -°°<t<  °° . 

b.  The  power  spectral  density  of  G(f)  exists. 
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a)   Analog  Signal 


-V; 


>   >  i   >   >   t  >   «   i   >   i   i   i 


■■— * 


•  i  h* 


b)  Digital  representation  of  analog  signal 


Figure  1.   Analog  and  Digital  Signals 
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c.   G(f)  equals  zero  for  frequencies  equal  to,  or 
greater  than  B.   B  is  the  highest  frequecy  which  can  be  re- 
solved and  f  is  defined  as  frequency. 

Further,  if  we  let  At  be  the  sampling  interval  in  seconds, 
so  that  At<l/(2B)  and  if  we  have  v(t)  sampled  at  intervals  of 
At,  giving  vi  samples,  i  =  -00,  .  .  ,0  ,1 ,  .  .  .  ,«°  ,  it  can  be  shown 
that  v(t)  can  be  reconstructed  uniquely.   So,  if  the  power  in 
v(t)  is  limited  to  a  band  less  than  V(2At)Hz,  then,  sampling 
at  an  interval  At  allows  v(t)  to  be  reconstructed  uniquely. 
Proof  of  this  augument  is  given  in  [Ref.   2].  In  this  way, 
the  requirement  of  sampling  at  least  twice   per  cycle  has  not 
only  been  established  but  is  entirely  sufficient. 
•  3 •   Limitations  on  Frequency  Resolution 

Due  to  real  world  limitations  of  not  being  able  to 
collect  a  signal  v(t)  of  infinite  length,  and  because  our 
actual  signals  are  not  necessarily  band-limited,  we  have  to 
modify  many  of  our  theoretical  assumptions.   Firstly,  we 
assume  we  can  get  a  long  record  length,  which  is  representative, 
at  least  over  the  range  of  frequencies  we  are  interested  in, 
of  the  signal  extending  in  time  to  plus  and  minus  infinity; 
secondly,  we  can  use  "low-pass"  filters  to  eliminate  un- 
wanted high  frequencies  before  the  signal  is  sampled. 

a.   Low  Frequency  Limitations 

If  the  signal  V(t)  contains  low  frequencies  it 
will  be  very  hard  to  distinqu'ish  then  in  the  interval  0<f<l/(2T), 
according  to  Rayleigh's  Criterion,  iHef.  21,  where  T  is  the 
record  length  in  seconds.   In  this  situation,  when  recording 
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a  finite  section  of  signal,  we  have,  in  fact,  truncated  the 
original  signal.   The  recorded  section  of  data  now  represents 
the  original  signal.   The  term  "block"  refers  to  a  further 
truncation  of  the  signal  into  smaller  sequential  sections  of 
data,  which,  when  added  sequentially,  will  give  a  truncated, 
sampled,  section  of  signal.   Thus  as  T  increases,  or  the 
length  of  signal  examined  increases,  the  lower  will  be  the 
frequency  which  can  be  resolved. 

b.  High  Frequency  Limitations 

The  highest  frequency  we  can  resolve  has  been 
established  as  one-half  the  sampling  rate.   In  other  words, 
at  least  two  samples  per  cycle  are  required.   The  high  fre- 
quency limit  is  more  commonly  known  as  the  "high  frequency 
cut-off  point,"  or  just  "cut-off  frequency,"  ( f c ) .   Sometimes, 
requirements  of  five  samples  per  cycle  are  set;  however,  this 
added  computational  problem  is  in  fact  unnecessary  and 
essentially  means  that  the  high  frequency  limit  of  analysis 
is  increased  by  a  factor  of  2.5.   If  1000  Samp/Sec  were  re- 
quired for  a  high  frequency  limit  of  500  cycles,  2500  Samp/ 
Sec  would  result  in  a  high  frequency  limit  of  1250  cycles. 
Usual  sampling  rates  vary  between  2  and  2.5  samples  per  cycle. 
In  this  study,  a  sampling  rate  of  two  samples  per  cycle  was 
used  for  the  turbulence  analysis.   This  rate  had  been 
established  as  satisfactory  by  Boston  (1970),  iRef.  lj . 

c.  Aliasing 

The  requirement  of  simply  sampling  the  highest 
frequency  of  interest  at  two  samples  per  cycle  is  adequate 
if  the  highest  frequency  of  interest  is,  in  fact,  the 
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highest  frequency  present  in  the  signal.   When  higher  fre- 
quencies are  present,  the  cut-off  frequency  becomes  the 
"folding"  or  "Nyquist"  frequency.   The  "folding"  comes  from 
the  fact  that  higher  frequency  energy,  "electrical  energy" 
in  our  case,  is  "folded"  back  into  lower  frequencies  *  around 
this  point.   To  eliminate  the  problem  it  is  necessary  to  either 
sample  at  higher  rates,  thus  moving  the  folding  frequency 
higher,  or  to  sharply  filter  the  signal  at  the  folding 
frequency . 

B.   FOURIER  TRANSFORMATION 

Many  studies  of  random  geophysical  processes  seek  to 

describe  the  distribution  of  energy  within  a  varying  process 

as  a  function  of  the  frequency  of  fluctuation  of  the  variable 

being  measured.   Sometimes  when  analysing  turbulence,  it  is 

desired  to  determine  the  distributation  of  energy  in  the 

rapidly  fluctuating  velocity  as  a  function  of  frequency.   The 

problem  is,  thus,  basically  one  of  transforming  data  from 

the  time-domain  into  the  frequency-domain. 

1 .   Fourier  Transformation  of  a  Continuous  Signal: 
Fourier  Integral 

One  of  the  most  often  used  techniques  for  this  trans- 
formation, from  the  time-domain  to  the  frequency-domain,  is 
through  the  use  of  the  Fourier  Intergral  Transform.   Accord- 
ing  to  this  procedure,  a  periodic  sinusoidal  signal  (shown  in 
analog  form  in  Figure  2a)  is  transformed  from  the  time-domain 
into  the  frequency-domain  through  the  transform-function 

V(f)=     /"v(t)e  "j27Tftdt  '  (1) 

where  v(t)  =  E  Sin  2tt  f0t  and  j=/^T 
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a)   Periodic  Sinosoidal  Signal 


b)  Fourier  coefficient 
of  transformed 
signal 
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c)   Complex  periodic  Signal 


d)  Fourier  coefficients 
of  transformed 
signal 


Figure  2.   Fourier  Transformations 
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This  transformation  is  shown  graphically  in  Figure  2b).   A 

complex,  almost  periodic  signal  composed  of  several  sine 

waves  of  differing  amplitudes  and  differing  frequencies  would 

be  similarly  transformed  into  the  frequency-domain,  as  shown 

in  Figure  2d) . 

2.   Fourier  Transform  of  Discrete  Data  Signal 

When  dealing  with  digitized  signals  or  discrete  data 

the  finite  form  of  the  Fourier-Transform  must  be  used: 

N-l    -j2TTfiAt 
V(fK)  .=  AtE   vi  e  _oo<f<oo  (2) 

where  v  is  a  complex  variable  and  K=0,...,N/2.   For  real 

data  vi  where  i=0,...,N-l,  the  sine  and  cosine  transforms 

becomes : 

N-l 
a(fK)  =  AtZ   v.  Sin  2iTf.At  (3a) 

i=0   x 
N-l 
b(f  )  =  AtZ    v.  Cos  27TfiAt  (3b) 

K      i=0 

Reference  2  suggests  that,  computationally,  this  would  follow 
a  flow  diagram  as  in  Figure  3.   The  value  of  f  is  computed  as 
a  function  of  the  number  of  data  points  and  the  length  of 
the  record  in  seconds  (T).   In  other  words,  increasing  the 
record  length  decreases  the  frequency  difference  between  co- 
efficients; or  a  longer  record  will  give  spectral  estimates 
closer  together  in  frequency.   Making  the  substitution  for 

f   the  sine  transformation  becomes: 
K 

N"1  IK 

a(KAf)  =  E   v±  Sin  2tt—  (4) 

i=0  N 

This  approach  to  computing  Fourier  coefficients  has  been 
limited  in  the  past,  due  to  the  computational  time  required. 
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Figure  3-   Flow  Chart  Showing  Fourier  Transform  Procedure 
for  use  with  Discrete  Data 
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The  number  of  calculations  required  increases  as  the  square 
of  the  number  of  data  points  increases.   The  integral  trans- 
form, in  the  past,  was  limited  to  studying  theoretical 
functions,  such  as  sine  waves  and  square  waves,  which  were 
relatively  well-behaved  functions.   For  cases  such  as  these, 
even  the  discrete  transform  could  be  used  to  readily  obtain 
the  Fourier-transform  of  a  signal.   With  the  introduction  of 
the  Fast  Fourier  Transform  (FFT),  algorithm,  computation  time 
has  been  greatly  reduced. 

3 .  Fourier  Transform  of  Truncated  Continuous  wave  Form 
If  a  signal  v(t)  exists  only  for  the  time  interval 

from  0  to  T  seconds,  and  is  zero  at  all  other  times,  its 

Fourier-transform  is: 

T     -j2irft 
V(f)  =  /  v(t)e      dt  (5) 

0 

If  v(t)  is  repeated  in  intervals  of  period  T  seconds,  the 
frequency  difference  between  coefficients  will  be  1/T  Hz.  The 
Kth  coefficient  will  be: 

T  -  j  2TT^t 

VK  =i  /  v(t)  e    T  dt  (oa) 

and  for  the  discrete  case:  ; 

N/2     -J2ir-|i 
V  =  L.Z   v*    e     T  Cob) 

K   T  i=1   1 

4 .  Convolution  of  Continuous  Signals 

The  finite-transform  of  a  finite-length  time  series 
can  be  viewed  as  the  product  of  a  finite-length  rectangular 
function  C    (t),  times  an  infinitely  long-time  history  v(t). 
The  finite  transform  of  v(t)  becomes: 

T  -j2TTft 

V(f)  =  /   CT/?(t)  v(t)  e       dt  (7) 

0   '  -  .  ■ 
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Since  products  transorm  into  convolutions,  the  convoluted 
Fourier-transform  becomes 

Vc(f)    =  V(f)x(JT/2(f)  (8) 

where 

T  .  -j2irft 

C|T/2   =   /   CT/2(t)    e  dt  (9) 

Figure  4(c)  shows  the  Fourier-transformation  of  the  rectangular 

function  and  Figure  4(d)  shows  the  convolution  of  a  sine  wave 

with  the  rectangular  function.   This  convolution  problem 

arises  when  switches  are  opened  and  closed  while  recording 

the  data.   The  side  lobes  of  the  convolved  function  can  be 

minimized  by  allowing  the  time  length  of  the  record  to  be 

sufficiently  large.   This  decreases  the  interval  0  to  1/T,. 

C.   POWER  SPECTRAL  DENSITY  FUNCTION 

1 .   Methods  of  Computing  Power  Spectral  Density 

There  are  three  methods  which  may  be  used  to  compute 
the  power  spectral  density  of  a  signal.   They  are: 

a.  Direct  Fourier  Transform  Method 

The  Fourier  transform  of  the  signal  is  computed 
and  from  this  the  mean  value  squared  is  determined. 

b.  Analog  or  Bandpass  Filter  Method 

The  signal  is  put  through  a  bank  of  bandpass 
filters  and  each  filter  output  is  squared  and  integrated. 

c.  Auto-correlation  Function  Method 

The  Auto-correlation  function  of  the  time  series 
is  computed,  and  then  its  Fourier-transform  is  computed. 
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a)   Rectangular  function 

(WW 


b)   Truncated  sine  function 


c)   Fourier  transform  of  rectangular  function 


^=7" 


d)   Fourier  transform  of  truncated  sine  function 


Figure  4.   Effect  of  Truncation  of  Sine  Function 
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Historically,  the  direct  method  was  developed  first, 
but  could  be  used  only  on  nicely  behayed  theoretical  signals. 
When  less  well  behaved  signals  where  analyzed,  smooth 
(deterministic)  theoretical  functions  were  replaced  by  discrete 
data  points  representing  the  signal.   This  necessitated  using 
the  discrete  Fourier-transform;  however,  for  large  data-sets, 
the  time  for  calculating  the  Fourier-transform  was  prohibitive. 
The  computer  program  FTOR  utilized  a  variation  of  the  direct 
method  employing  the  Fast  Fourier  Transform  (FFT). 

Using  the  sampling  rates  1000-5000  samp/sec,  (SPS), 
and  having  the  ability  to  select  the  block  sample  length  and 
the  number  of  blocks  desired  for  a  particular  run,  no  pro- 
blems were  encountered  in  which  filler  data,  usually  zeros, 
had  to  be  inserted  into  a.  block.   The  block  size  was  always 
chosen  as  an  integral  power  of  two. 

2 .   Typical  Power  Spectral  Density  Functions 

The  Power  Spectral  Density (PSD)   plot  for  random  data 
shows  the  distribution  of  electrical  power  within  the  signal 
as  a  function  of  frequency.   Several  characteristic  PSD  plots 
are  encountered.   Figure  5(a)  is  a  PSD  plot  of  a  pure  sine 
wave.   It  gives  the  Dirac-delta  function  which  implies  that 
the  power  at  the  sine  frequency  is  infinitely  large,  and  zero 
at  all  other  frequencies.   Figure  5(a)  is  a  PSD  of  a  Gaussian' 
random  signal.   The  PSD  of  this  signal  is  constant.   Figure 
5(c)  shows  the  PSD  of  a  random'  signal  carrying  a  sine  function 
on  it.   The  PSD  plot  is  the  sum  of  the  power  spectra  of  the 
random  signal  and  sine  figured  separately. 
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c)   Gaussian  Signal  Plus  Sine  Wave 


Figure  5.   Characteristic  PSD  Plots 
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3.  Problem  of  Single-Sided  Spectrum 

The  Fourier  coefficients  under  the  transformation 

-J2itft' 
Vlf.)  =/_«>  v(t)e      dt  (10) 

exist  for  both  positive  and  negative  frequencies.   The  co- 
efficients can  be  transferred  to  a  single-sided  frequency  plot 
by  simply  doubling  each  coefficient  as  it  is  plotted  only  a- 
long  the  positive  frequency  axis.  If  S(f)  represents  the  two- 
sided  spectral  density  function  and  G(f)  represents  the  single- 
sided  spectral  density  function,  the  single-sided  spectral 
density  function  equals  twice  the  double-sided  density  function 
G(f)=2S(f).   Likewise,  if  a  watt  meter  was  used  to  find  the 
power  in  a  signal  as  a  function  of  frequency,  the  values  ob- 
tained would  have  to  be  divided  by  two  before  plotting  on  a 
two-sided  spectrum. 

4 .  Power  and  Energy  Signals 

The  average  ■  electrical  power  in  a  fluctuating 
voltage  signal  v(t)  is  given  by 

P  =  lim  _l/a  |v(t)|2  dt  (11) 

a-*-*  2  a 

Mix  [Ref.  3]  defined  a  power  signal  as  one  which  is,  for  all 
practical  purposes,  infinitely  continuous.   Energy  signals 
are  pulse-like  in  form  and  are  given  by: 


E  =  /"  |v(t)|2  dt  .  (12) 

—  00 
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The  power  in  a  periodic  signal  is  shown  from 
Parseval's  Theorem  to  be: 

P  =  i^;+T  ivct)i2  dt  (13) 

Thus,  the  power  in  a  continuous  sine  function,  v(t)=sin  377t, 
would  be 

P  =_1  /T  Sin2  377t  dt  =  1/2  watt  (14) 

T   0 

If  a  real  wide  band  signal  were  to  be  analyzed  by  using  a  tun- 
able filter  of  band  width  of  6  f ,  a  plot  of  power  versus 
frequency  as  shown  in  Figure  6  might  result.   Here,  the  real 
power  from  zero  to  an  upper  frequency  was  determined,  divided 
by  two  and  the  values  folded  back  into  negative  frequencies. 
5.   Power  from  the  Fourier  Coefficients 

Parseval's  Theorem  leads  to  the  definition  of  power 
in  terms  of  the  Fourier  coefficients 


p  =  ^/*:+T  Mt)i2  dt  =  i  jvKiz  d5) 


Thus,  knowing  v(t),  the  average  power  can  be  computed.   The 
average  power  in  a  one  volt,  60Hz  sine  wave  is: 

v(t)  =  Sin2  60t 

P  =  i/^  Sin2  377t  dt  =  1/2  watts  (16) 

TO 

Using  the  Fourier  coefficients: 

V.-l  =  1/2   V-l  =  1/2  VK>1  =  0  . 

P   =   Z  |VJ2    =   1/4    +   1/4    =    1/2  (17) 

K=-«>        K 
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Figure  6.   Two -Sided  Power  Spectrum  of  Wide  Band  Signal 
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Figure  7.   Power  Spectrum  of  1  volt  60  Hz  Sine  Signal 
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If  a  watt  meter  was  used  to  check  the  power  of  this 
sine  signal,  one-half  watt  would  be  found  at  60Hz,  and  zero 
watts  at  other  frequencies.  If  a  plot  of  power  versus  fre- 
quency were  made,  one  similar  to  Figure  7  would  result.  Since 
the  values  have  been  plotted  for  positive  and  negative  fre- 
quencies, the  total  power  has  been  reduced  by  one-half  and 
the  values  have  been  plotted. 

D.   FAST  FOURIER  TRANSFORM 

1 .   Computational  Economy  Afforded,  by  FFT 

The  Fast  Fourier  Transform  (FFT)  algorithm  is  a 
fast  method  for  computing  the  finite  Fourier  transform  of  a 
series  of  N  data  points.   The  FFT  computation  requires  N 

log?  N  computations,  which  leads  to  a  substantial  saving  in 

2 
computer  time  over  the  old  method  which  required  N  operations. 

The  FFT  economy  of  computational  effort  is  greatest  for  large 
values  of  N. 

Reference   [4]  gives  the  general  computational 
routine  for  figuring  the  FFT.   It  is  pointed  out  that  the 
FFT  is  general  in  that  N  is  not  necessarily  a  power  of  two; 
however,  by  selecting  N  to  be  a  power  of  two,  further  computa- 
tional savings  result.   When  this  requirement  is  met,  the 
FFT  algorithm  is  essentially  a  successive  doubling  operation. 
Thus  for  N=2J ,  only  Nj  multiply-add  operations  would  be  re- 
quired under  the  FFT  assuming  the  necessary  complex  exponential 
table  of  values  has  been  computed  in  advance. 
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2.   Importance  of  FFT  in  Turbulence  Analysis 

Computationally,  the  FFT  is  an  indispensable  aid  in 
the  PSD  investigation  of  turbulence.   Due  to  the  extremely 
high  sampling  rates  necessary  for  studying  high  wave  number 
processes,  enormous  amounts  of  data  are  collected.   A  five 
minute  section  of  a  fluctuating  temperature  signal  sampled 
at  a  rate  of  4000  samples/second  (SPS)  generated  almost  a 
million  and  a  half  data  samples.   Using  the  old  method  for 
computing  the  Fourier  transform  for  this  quantity  of  data, 
several  thousand  billion  operations  would  be  required.   Clearly 
the  time  element  in  this  procedure  would  prohibit  the 
calculation  on  any  but  the  fastest  computers.   Turning  to  the 
FFT  the  operation  would  only  take  about  21  million  operations 
or  an  improvement  in  excess  of  5  orders  of  magnitude  in 
computational  time.   The  FFT  thus  brings  the  study  of 
turbulence  signals  within  the  bounds  of  computational 
feasibility . 

The  ability  to  determine  the  power  spectral  density 
function  which  is  computed  from  the  Fourier-transform  and 
which  is  such  an  important  aspect  of  turbulence  analysis  is 
made  possible  through  the  use  of  the  FFT  algorithm. 
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III.  NAVAL  POSTGRADUATE  SCHOOL  DIGITIZATION  FACILITY  AND 

SPECTRAL  ANALYSIS  PROGRAMS 

A.   NPS  COMPUTER  FACILITIES  FOR  DIGITIZATION  AND  PSD  ANALYSIS 

The  Naval  Postgraduate  School  has  two  sophisticated  computer 
systems  which  are  used  in  the  digitization  and  PSD  analysis 
procedures.   One  is  a  Hybrid  system  which  consists  of  an 
analog  computer,  COMCOR  Ci  5000,  which  is  electrically  inter- 
faced with  a  digital  computer,  XDS  9300.   The  XDS  9300  con- 
tains 3^K  bytes  of  core  storage  and  uses  an  octal  number  base. 
This  base  consists  of  binary  numbers  made  up  of  3-bit  digits. 
The  second  system  is  an  IBM  360/67  which  has  a  core  storage  of 
762K  bytes.   It  uses  a  hexidecimal  number  base  which  consists 
of  binary  numbers  of  4-bit  digits. 

1 .  Hybrid  Computer 

The  Hybrid  computer  facility  is  located  on  the  fifth 
floor  of  Spangel  Hall.   The  equipment  disposition  is  shown 
in  Figure  8.   Though  the  individual  user  performs  all  computer 
operations,  technical  assistance  is  availabe  from  the 
Electrical  Engineering  Computer  Laboratory  staff.   This 
facility  is  used  for  the  analog-to-digital-  conversion  of 
signals. 

.......  * 

2.  IBM  360/67  Computer 

The  large,  digital  computation  facility  is  located  on 
the  ground  floor  of  Ingersol  Hall.  Computer  Center  personnel 
handle  all  computer  operations.   Individual  programs  are  input 
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to  the  computer  through  a  card  reader  and  requests  for  the 
subsequent  input  of  programs,  tapes,  disk  memory,  etc.  are 
submitted  "over  the  counter"  to  computer  personnel.   A  24- 
hour  operating  schedule  is  usually  maintained.   All  PSD 
analyses  of  digital  tapes  were  performed  on  this  computer. 

B.   HYBRID  COMPUTER  SYSTEM 

Computer  time  on  the  Hybrid  computer  system  is  signed  for 
in  advance  in  the  Electrical  Engineering  Computer  Laboratory. 
The  whole  system  is  operated  on  a  self-service  basis;  however, 
the  computer  center  staff  is  available  to  assist  during 
working  hours.  ■  Computer  time  was  easiest  to  obtain  early 
in  the  term  or  between  terms  when  the  work-load  was  lightest . 
The  facility  is  available  24  hours  a  day  and  as  one  becomes 
more  proficient  in  equipment  operations,  night  or  week-end 
operations  can  be  scheduled  if  week-day  operations  are  precluded 

1.   The  Analog  Computer  (Ci  5000) 

The  Analog  Computer  contains  the  input  points  for 
raw  signals  input.   It  also  has  control  features  for  signal 
amplification,  selection  of  sampling  rate  and  starting  and 
stopping  of  the  digitizing  procedure.   The  two  removable 
patch  boards  are  the  heart  of  the  system.   Individual 
patchboards  have  been  reserved  solely  for  analog-to-digital 
conversion  (board  #24). 
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The  logic  board,  the  smaller  of  the  two  patchboards, 

was  utilized  to,  set  the  sampling  rate  of  the  analog  computer. 

The  logic  board  with  its  patch  is  shown  in  Figures  9  and  17. 

The  sampling  rate  was  changed  by  turning  a  counter-like  number 

display,  below  and  to  the  left  of  the  logic  board,  as  shown 

in  Figure  17.   The  counter  operated  a  voltage  divider  which 

was  electrically  connected  to  a  resistance  input  on  the  logic 

board.   The  resistance  value  at  this  point  was  changed  by 

inserting  a  resister  in  one  of  4  positions.   The  resistance 

value  was  then  divided  by  the  counter  setting  plus  one  to 

give  the  sampling  rate  (in  samples  per  second)  desired. 

Example:   lOOkc   =  4Kc  or  4000  samples  per  second. 
(24+1) 
No  external  leads  are  connected  to  the  logic  board. 

. The  analog  patchboard  is  the  input  point  for  all 

external  signals  entering  the  analog  computer.   The  input 

signal  from  a  tape  recorder  or  signal  generator,  after  passing 

through  signal  conditioning  equipment  is  input  to  one  of  the 

analog  board  amplifiers.   (Figure  15  shows  inputs  going  into 

amplifier  A001).   Various  gain  factors  can  be  applied  to  the 

input  signal  to  bring  it  up  to  a  value  optimum  for  utilization 

of  the  analog  computers  dynamic  range.   Figure  10  shows  the 

keyboard  of  the  analog  computer  which  is  used  to  control  the 

desired  operational  mode  of  the  analog  computer. 

.    ...  x 

2.   Digital  Computer  QCSD  9300) 

The  heart  of  the  digital  computer  system  is  the  Xerox 
Data  System  9300  Central  Processor  Unit  (CPU).  Two  tape  drive 
units,  line  printer,  card  reader  and  teletype  unit  are  inter- 
faced with  the  CPU.   Figure  11  shows,  diagrammatically ,  the 
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equipment  used  in  digitization.   The  tape  driye  units  are 
used  for  data  input  from  tape  to  computer  or  for  data  out- 
put from  computer  to  tape;  the  line  printer  can  be  used  for 
program  listings  and  printing  results;  the  card  reader  as  pro- 
gram input  for  compiling  programs;  the  teletype  unit  for  input 
of  parameters  necessary  to  control  the  "multi-channel  analog- 
to-digital  Program."  The  CONTROL  CONSOLE  of  the  XDS  9300  is 
used  to  compile  and  to  initialize  the  program.   After  the  pro- 
gram has  been  compiled,  the  teletype  is  used  to  select  one  of 
seven  digitizing  options.  The  digital  and  analog  computers  are 
connected  through  the  Analog-to-Digital  Converter  (ADC),  (Fig. 12). 
3 .   Multi-Channel  Analog-to-Digital  Program 

Several  programs  have  been  developed  by  the  computer 
laboratory  staff  to  control  the  XDS  9300  system  during  multi- 
channel digitizing  operations.   A  card  deck  of  the  latest  re- 
vision of  the  program  is  maintained  in  the  computer  lab.  This 
program  can  be  input  from  either  cards  or  from  a  tape  which 
has  had  the  machine  language  program  stored  on  it .  The  com- 
piling sequence  takes  about  five  minutes  using  the  card  input 
and  only  about  30  seconds  using  the  tape  input. 

Once  the  program  has  been  compiled,  one  of  the  seven 
Program  Control  Options  can  be  selected  by  typing  one  of  the  , 
options  followed  by  carriage  RETURN  (C/R)  on  the  teletype 
(see  Figure  13.)-   The  Program.  Control  Options  are: 

1.  Enter  new  parameters.  (NSAMP,  NCHAN,  NREC,  ITAPE) 

2.  Start  digitizing  the  analog  input  signals 

(Digitization  actually  starts  when  manual  switch  DSI  on  Ci  5000 

is  thrown  to  up  position) . 
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Figure  12.   Block  Diagram  of  Analog-to-Digital  Conversion 
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Figure    13-      Teletype 
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3.  Write  End-of-File  on  tape. 

4.  Rewind  tape  to  load  point. 

5.  Skip  files  (the  number  of  files  to  be  skipped 
specified  by  teletype  input  of  a  four  digit  number,  i.e.  skip 
five  files  000  5  as  shown  in  Figure  14). 

6.  Print  digitized  data  (the  number  of  lines  to  be 
printed  specified  by  teletype  input  of  a  four  digit  number, 
i.e.  print  ten  lines  of  digitized  data  0010. 

7.  Actuate  the  Digital-to-Analog  subroutine  for  the 
next  single  block  of  data  encountered  and  input  the  data  into 
the  strip  chart  recorder  (a  subsequent  computer  comment  types 
"start  strip  recorder  and  type  C/R  -  carriage  return").   The 
strip  recorder  must  first  be  connected  to  the  digital-to-analog 
output  on  the  analog  patchboard  on  the  Ci  5000. 

After  the  multi-channel  digitization  program  has  been 
compiled  and  the  input  light  on  the  teletype  is  on,  tape 
parameters  input  through  the  teletype  are: 

NREC  -  '  sets  the  maximum  number  of  records  to  be 
digitized  (see  Figure  14). 
'  NSAMP  -   sets  the  number  of  digitized  samples  contained 
in  each  record  on  the  seven-track  tape 
( see  Figure  14 ) . 
ITAPE  -   sets  the  magnetic  tape  unit  into  which  the  ' 
digitized  samples  are  sent.   This  number  must 
be  the  same  as  the  dial  number  selected  on 
the  front  of  the  tape  drive  unit  (see  Figure  14) 
NCHAN  -   sets  the  number  of  analog  channels  to  be 
digitized  (see  Figure  14). 

42 


a} 

-p 

•H 
hO 

•H 

Q 

I 

o 
-p 

I 

hO 
O 
i— I 

CC 

C 

< 

O 

<H 

CO 
P 

a 


0)  i-i 
D,  O 

-P     P 

CD  C 
rH  O 
CD  O 

£ 
ctf 
Ih 
50 

o 


ctf 
a 

•H 

& 

Eh 


CD 

•H 
Ph 


^3 


C.   ANALOG-  TO-DIGITAL  OPERATIONS 
1 .   Equipment  Set  Up 

a.  Analog  Tape  Recorder 

The  most  commonly  used  record/playback  device  in 
the  Oceanography  Department  is  the  Sangamo  tape  recorder 
Model  3562,  which  can  record  up  to  14  tracks  of  data,  and 
two  edge-tracks  of  voice  on  one  inch  magnetic  tape.   Two  modes 
of  operation  are  possible:   direct  and  FM.   It  is  general 
practice  to  only  use  the  FM  mode,  which  records  from  D.C.  to 
some  upper  limit  depending  on  the  tape  speed.   The  FM  mode  is 
most  noise  free.   Direct  electronics  are  used  for  higher 
frequency  data  when  no  D.C.  information  is  required. 

The  magnectic  heads  and  tape  rollers  of  the  tape 
recorder  should  be  cleaned  with  isopropyl  alcohol  after  pro- 
longed use  to  reduce  noise  during  playback. 

b.  Filters 

KHRON-HITE,  Model  3321  filters  have  been  most 
widely  used  in  the  past  for  filtering  out  unwanted  high  and 
low  frequencies  from  a  signal.   They  can  either  by  connected 
singlely  in  a  low-pass  mode  or  in  series  as  a  band  pass  fil- 
ter.  After  the  filter  mode  has  been  selected  the  raw  signal 
is  input  to  the  filter  and  the  filter  output  connected  to  the 
analog  patchboard  as  an  input  to  the  analog  computer. 

c.  Analog  Patchboard 

The  analog  board  should  be  patched  as  shown  in 
Figure  15  if  two  channels  are  to  be  digitized.   If  only  one 
channel  is  desired  the  signal  input  should  be  input  into 
amplifier  A001  only.   Gain  factors  may  be  patched  as  desired. 
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The  outputs  of  each  amplifier  go  to  both  the 
oscilloscope  input  points  and  the  ADC  input  points  (T500  and 
T501).   If  an  analog  playback  of  the  digitized  signal  is 
desired,  inputs  from  T420  and  T^21  (see  Figure  15)  must  be 
made  into  the  recorder  inputs  (1-8). 

d.  Logic  Patchboard 

The  patch  for  the  logic  board  is  shown  in  Figure 
17.   This  is  for  continuous  mode  digitizing  in  which  digitization 
is  started  by  throwing  switch  DS1 .   The  only  logic  setting 
necessary  is  the  selection  of  the  required  resistance  value 
(Figure  17)  which  in  conjuction  with  the  wheel  counter  sets 
the  sampling  rate. 

e.  Oscilloscope 

The  "scale  illumination"  dial  energizes  the  oscil- 
loscope (see  Figure  16).   The  sampling  pulse  was  usually 
patched  into  "channel  1"  by  throwing  switch  Yl  to  the  left. 
This  permitted  continuous  monitering  of  the  sample  pulse  fre- 
quency and  width.   Dial  DFOO  (see  Figure  17)  should  be  set  to 
Ms.  1  and  the  width  of  the  sample  pulse  adjusted  with  this  dial. 
2 . t     Energizing  the  Ci  5000  Computer 

If  the  power  on  the  CI  5000  has  been  secured  (see 
power  light  in  Figure  9b)  it  can  be  turned  on  by  depressing 
the  "on"  switch.   Next,  the  buttons  KEYBOARD,  POTSET  and 
RESET  are  depressed  in  that  order  (see  Figure  10). 
3.   Energizing  the  JSDS  9300 

The  step  by  step  procedure  below  should  be  followed 
in  energizing  the  Hybrid  computer  system.   The  latest  update 
to  this  procedure  is  kept  on  the  XDS  9300  Control  Console. 
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Figure  17.   Ci  5000  Logic  Board  and  Operating  Switches 
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a.  XDS  9300  Power  on  (if  power  is  off),  see  Figure  18 

1)  Depress  IDLE 

2)  With  RESET  depressed  press  POWER 

3)  Turn  teletype  switch  to  ON 

4)  Press  CLEAR  and  CLEAR  FLAGS  simultaneously 

b.  Energizing  the  line  Printer  (if  the  READY  light  isn't  on) 

1)  Insure  that  POWER  button  on  line  printer  is  on. 
If  not,  press  POWER  button. 

2)  When  lower  half  of  POWER  is  lit  (within  one  min. 
after  POWER  pressed)  press  READY,  NOTE:  READY  must  be  off  to  ad- 
vance paper.  A  whole  page  is  advanced  with  TOP  of  FORM  and  the 
paper  is  advanced  a  single  line  with  SINGLE  SPACE. 

c.  Card  Reader  (if  A/D  program  being  input  from  cards) 

1)  Place  BOOT  card  in  front  of  the  JOB  card  in  the 
A/D  deck. 

2)  Put  cards  into  card  reader  (see  Figure  19)  face 
down,  top  outward  (toward  you)  so  that  column  1  is  first  into  reader. 

3)  Put  card  reader  press  in  order:   POWER  and 
START . 

d..  Program  Compile  and  Load 

1)  Ready  the  Line  printer  (sequence  b  above). 

2)  With  cards  in  place  ready  the  card  reader 
(sequence  c  above). 

3)  At  9300  console  turn  off  (by  depressing)  all 
SENSE  switches  which  are  lit . 

4)  Press  in  order:   IDLE;  RESET:  RUN:  CARDS. 

5)  A  JOB  is  printed  out  by  the  teletype  indicat- 
ing compiling  has  begun.  "' 
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6)   Upon  compilation,  the  teletype  will  type 
END  OF  ASSEMBLY.   Shortly,  thereafter,  the  INPUT  light  will 
go  on,  indicating  the  tape  parameters  can  be  typed  in. 

e.  Equipment  NOT-READY  Condition 

1)  If  a  device  is  not  ready  when  needed,  it  will 
be  called  out  by  a  teletype  message.   Simply  making  the  device 
ready  will  clear  the  pause  in  the  operation. 

2)  In  the  event  of  a  FEED-CHECK  failure,  inspect 

the  card  that  caused  the  problem  (the  one  on  the  bottom  of  the 
input  hopper).  .  If  the  card  is  found  to  be  damaged  on  the 

leading  edge,  repunch  the  card  and  replace  it.   Put  the  un- 
read cards  back  in  the  hopper.   Ready  the  card  reader  by  press- 
ing RESET  and  START. 

3)  Line  printer  not-ready  may  be  caused  by  its 
being  out  of  paper. 

f.  XDS  9300  POWER  OFF  (after  normal  working-hours 
only  ) 

1)  Press  IDLE. 

2)  Turn  teletype  switch  to  OFF. 

3)  Press  IDLE. 

4)  With  RESET  depressed,  press  POWER. 
4 .   Mounting  Magnetic  Tapes 

It  is  very  important  to  use  currently  certified  tapes 
which  are  free  of  permanent  errors.   Scratches,  tears,  oil 
spots,  etc.  cause  tape-writing  errors.  The  only  recourse  is  to 
write  an  END  OF  FILE  on  the  tape  at  the  end  of  the  bad  file 
and  commence  digitizing  on  the  next  section  of  tape. 
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SCIENTIFIC  O&TA  SYSTEMS 


Figure  20 


7-Track  Magnetic  Tape  Drive  Unit 
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Figure  2  0  shows  one  of  the  two  tape  driye  units  which 
may  be  selected  for  tape  mounting.  Turn  the  mode  dial  to 
MANUAL  READ.   Before  mounting  the  tape  a  circular  plastic  file 

WRITE  RING  must  be  inserted  into  the  circular  slot  on  the  back 
of  the  seven-track  tape  being  mounted.   Mount  the  tape.   If 
the  ring  is  in  place  the  FILE  PROTECT^ Figure  21  light  on  the 
top  row  of  indicator  lights  will  go  out,  indicating  that  data 
can  be  written  on  the  tape.   Tape  threading  instructions  are 
posted  on  the  inside  of  the  protective  doors  enclosing  the 
tape  reels.   Once  the  tape  is  secure  around  the  take-up  reel, 
it  should  be  wound  around  at  least  six  more  times  to  prevent 
the  tape  from  slipping  off  the  take-up  reel,  when  the  main  se- 
curing latch  is  closed.   With  the  securing  latch  in  place  and 
the  tape-idler-arm  down,  depress  FORWARD-DRIVE.   The  tape 
will  move  forward  and  stop  when  it  reaches  the  metallic  strip. 
called  the  LOAD  POINT.   The  indicator  light  LOAD  POINT  on  the 
top  row  of  indicator  lights  will  go  on.   After  the  tape  has 
been  mounted  and  advanced  to  the  LOAD  POINT,  the  tape  density 
is  selected  by  turning  the  DENSITY  dial  to  either  256  or  556. 
This  determines  the  recording  density  in  bits  per  inch  and  is 
normally  set  to  556.   The  TAPE  UNIT  dial  is  turned  to  any  one 
of  several  numbers,  i.e.  1.   (The  teletype  input  parameter 
ITAPE  must  match  this  number,  i.e.  ITAPE  =1.)   Finally,  the 
mode  dial  is  turned  to  automatic.   The  TAPE  DRIVE  UNIT  is 
ready. 
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5.   Variable  Tape  Digitizing  Parameters 

Once  the  energizing  and  sampling  sequences  have  been 
completed  (the  INPUT  light  on  the  teletype  should  be  on),  the 
tape  options  are  entered  through  the  teletype  by  typing  a 
single  digit  number  and  then  depressing  the  CARRIAGE  RETURN., 
C/R  button.   In  order  to  input  new  parameters,  "1  C/R" 
would  be  typed. 

Then  the  following  example  parameters  could  be  typed 
in: 

NSAMP  =  2048  (C/R) 

NREC   =  100   (C/R) 

NCHAN  =  2,  ITAPE  =  1,  NDEL  =  2000  (C/R) 
The  computer  would  then  type  OPTION  (II),  soliciting  a  further 
response  from  the  user.   If  analog-to-digital  collection  was 
to  commence,  the  programmer  would  type  "2  (C/R)."   If  the 
Analog  computer  switch  DS0  was  in  the  UP  position,  digitization 
would  begin  when  DSI  was  thrown  to  the  UP  position. 

As  would  be  expected,  a  high  sampling  rate  would  fill 
up  100  records  faster  than  a  lower  sampling  rate.  Thus,  if 
a  sampling  rate  of  2000  SPS  had  been  selected,  and  the  above 
tape  parameters  had  been  selected,  the  digitization  would 
result  in  204, 800  samples  being  collected  and  a  total  signal 
length  of  102.4  seconds. 

2048  SAMPLES  x    ioo  RECORDS  =  204,800  SAMPLES 

204,800  SAMPLES  ,  e-1  ?  nr?nrmuc, 

2000  SAMPLES/SEC.  x  2  CHAN    0±  ' d    ^ouMJb 


5'6 


There  would  only  be  102,400  samples  of  each  channel;  however, 
there  were  two  channels  being  digitized  simultaneously.   Figure 
22  shows  how  digitized  data  for  one  and  two  channel  digitization 
would  be  formatted  on  a  seven-track  tape. 

D.   CONVERT  PROGRAM 

The  CONVERT  PROGRAM  converts  the  seven-track  octal  base 
data  samples  into  nine-track  hexidecimal  samples.   The  basic 
character  of  the  data  samples  being  in  blocks  on  the  magnetic 
tape  is  maintained;  however,  two  additional  parameters   (see 
Figure  23)  are  affixed  to  the  beginning  of  each  block  by  the 
CONVERT  program.   These  are  the  numerical  values  KMAX  which 
is  set  equal  to  the  maximum  number  of  samples  per  block  and 
NCHAN,  the  number  of  channels  of  analog  data  digitized  on  the 
seven-track  tape.   This  change  in  the  block  format,  the  change 
in  number  base  and  the  addition  of  two  new  parameters  is 
necessary  to  put  the  tape  in  the  proper  format  for  input  into 
the  FTOR  program. 

Figure  24  shows  the  CONVERT  procedure  in  flow-chart  form. 
The  program  processes  all  of  the  records  in  a  file  and  the 
number  of  files  to  be  converted  is  specified  in  the  IF  state- 
ment which  compares  the  number  of  files  completed  with  the 
number  specified  for  processing.   The  number  in  parenthesis 
is  set  to  the  number  of  files  to  be  converted.   The  JOB  CONTROL 
LANGUAGE  (JCL)  cards  which  follow  the  CONVERT  PROGRAM  specify 
which  files  on  a  multi-file  tape  are  to  be  processed.   Jones 
[Ref.  5  Appendix  D]  gives  the  typical  JCL  cards  used  to  con- 
vert five  files. 
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Figure  22. 


7-Track  Tape  Formatting  for  Single  and  Dual 
Channel  Digitization 
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Figure  23.   9-Track  Tape  Formatting  for  FTOR  Program 
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Figure  2k.      CONVERT  Flow  Chart 
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E.   NAVAL  POSTGRADUATE  SCHOOL  PSD  COMPUTER  PROGRAMS 

Several  computer  subroutines  are  -available  on  the  IBM 
360/67  for  computing  the  Fourier  transform  of  digitized  in- 
put signals.   Most  programs,  however,  are  designed  for  the 
input  of  relatively  small  data  sets,  usually  from  Hollerith 
cards.   These  programs  could  be  designed  to  compute  the 
Fourier  transform  of  large  amounts  of  data  on  magnetic  tapes; 
however,  a  very  powerful  series  of  programs  is  available  in 
the  Naval  Postgraduate  School  computer  library  which  will 
compute  the  FFTof  large  data  sets.   These  programs  are  UBCFTOR 
UBCSCOR  and  UBCFCPL  which  are  stored  under  files  one,  two  and 
three  respectively  on  tape  NPS  216.   The  only  particular  format 
requirement  for  input  data  is  that  the  digitized  samples  must 
be  arranged  in  groups  of  an  integral  power  of  two  samples  in 
each  "record."   A  record  is  2$    samples  where  j  is  an  integer.. 
Many  records  may  be  found  on  a  single  magnetic  tape  (see 
Figure  22).   Since  this  FFT  program  package  is  especially 
well  suited  for  analyzing  large  digitized  data  sets  stored  on 
magnetic  tape,  it  serves  as  an  essential  component  of  the 
Naval  Postgraduate  School's  signal  processing  capability. 
Figure  25  shows  schematically  the  complete  digitization  and 
PSD  analysis  procedure  using  these  programs. 

The  following  description  of  the  FFT  package  was  taken 
from  Dobson  [Ref.  6] .   Editorial  changes  were  made  to  update 
the  information  with  the  Naval  Postgraduate  School  facility. 
The  programs  were  written  by  J.  F.  Garett  and  J.  R.  Wilson 
while  students  at  the  Institute  of  Oceanography  at  the 
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Figure  25. 


Schematic   Diagram  of  Complete  Digitization 
and  PSD  Analysis  Procedure 
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University  of  British  Columbia.   The  UBCFTOR  program, 
written  in  Fortran  IV,  instructs  the  IBM  360  to  read  the 
digital  data  tapes,  store  a  block  of  data  into  the  computer 
memory,  compute  the  FFT  of  the  data,  and  then  store  the 
resulting  Fourier  coefficients  on  another  tape.   From  these 
coefficients  the  UBCSCOR  program  computes  the  spectral  density 
for  single  channel,  or  the  cross-spectral  density  between 
selected  pairs  of  channels.   The  spectral  density  values  may 
be  "output"  in  either  tabular  form  or  in  graphic  form  on  the 
Calcomp  plotter.   The  UBCFCPLT  plots  the  amplitudes  of  the 
Fourier  coefficients  as  a  function  of  frequency. 
1.   Fourier  Coefficient  Program  UBCFTOR 

This  program  calls  several  subroutines  which  read 
in  data  from  the  digital  tapes  generated  in  the  A/D  phase. 
The  maximum  number  of  data  points  which  can  be  stored  in  the 
computer  memory  under  this  program  are  8192  or  2    ->    samples. 
This  restriction  was  imposed  by  PKFORT,  the  subroutine  which 
computes  the  FFT .   The  details  of  this  subroutine  are  listed 
under  PKFORT  in  the  computer  library.   For  optimal  efficiency 
PKFORT  has  been  designed  to  read  in  data  blocks  which  con- 
tain an  integral  power  of  two  (i.e.  2^  )  samples.   PKFORT  is 
then  called  and  the  Fourier  transform  is  computed,  which  re- 
sults in  2^~   Fourier  coefficients.   Since  only  one  block  of 
data  is  analyzed  at  a  time,  the  original  signal  is  divided 
into  short,  sequential  sections  of  signal  and  a  new  length 
of  record  T  transpires.   The  resulting  coefficients  are  stored 
on  another  magnetic  tape  and  they  become  the  basis  for  the 
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spectra.   For  this  block  of  data,  the  highest  frequency 
for  which  a  coefficient  has  been  computed  is  the  Nyquist  fre- 
quency fs/2,  where  fs  is  the  sampling  frequency.   The  lowest 
frequency  which  can  be  resolyed  for  this  data  block  is  1/T 
or  2^/fs. 

The  Fourier  coefficients  are  written  into  blocks  on  a 
"coefficient"  tape  which  also  contains  identifying  information 
such  as  the  block  number,  sampling  frequency,  number  of  samples 
in  each  block,  etc.   Upon  completion  of  the  transformation  and 
writing  sequence  for  one  block,  the  next  block  of  data  is  read 
from  the  digitized  tape.   The  sequence  is  repeated  for  as  many 
blocks  as  specified  on  the  input  data  cards,  until  an  end-of- 
file  mark  is  sensed  on  the  digitized  tape  or  until  a  blank 
card  is . encountered  in  the  input  data  cards. 
2.   Spectral  Analysis  Program  -  SCOR' 

Once  the  Fourier  coefficients  have  been  computed, 
the  next  step  is  to  convert  these  values  into  spectral  values. 
The  program  SCOR  reads  the  coefficient  tape  generated  by  FTOR 
and  averages  the  PSD  values  over  32  bandwidths..   The  data 
card  section  of  the  SCOR  program  specifies  the  number  of 
channels,  the  number  of  blocks  to  be  analyzed  within  a  file, 
the  block  number  where  the  analysis  is  to  begin,  whether  the 
spectra  for  a  single  channel  or  the  cross-spectra  between 
channels  is  to  be  computed,  the  type  of  bandwidth  desired 
(constant  or  logarithmic),  etc.   Just  as  with  FTOR  procedures, 
SCOR  reads  and  analyzes  only  one  block  of  coefficients  at  a 
time . 
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If  smoothing  of  the  coefficients  is  desired,  various 
smoothing  functions  may  be  selected.   The  "harming"  option 
performs  a  three  point  running  average  on  the  data  with  weights 
of  -1/4,  1/2,  -1/4.   References  2  and  6,  contain  further  de- 
tails of  the  smoothing  functions. 

If  individual  Fourier  coefficients  of  two  channels  are 

1/2 
R-,  +  il,  and  R  +  ilp,  where  i=  (-1)    and  x  equals  the  block 

length  in  seconds  (  and  6F=l/x  where  <Sf  is  the  bandwidth  between 

Fourier  coefficients)  then  the  power  spectral  density  is  given 

by: 

x(Rf  +  if)      Rf  +  if 

*11(f)=  ± i-  =  .-i ± 

-  2  26f 

/r,2    T2N       2     2 

,        4   t(R0  +  I0)      r0  +  I0 

and   4,   (f)=  _^ L_  =   ^2 2 

2  25f 

The  cross-spectral  values  are  given  by  the  co-spectrum: 

C  (f)  =  t(R1R2  +  I1I2)  _  R1R2  +  ^2 
0  2  2<Sf 


and  the  quad-spectrum: 

t(r  r   _  it  )     R  R  -  R  R 

fs»  2  1     2  1       2  1     2  1 

Qu1P(f)  =  =  

ld  2  26f 


The  factor  of  2  in  the  above  equations  makes  the  integral 
under  the  power  spectrum  over  positive  frequencies  equal  to 
the  signal  variance.   Phase  corrections  can  be  made  at  this 
point  to  correct  for  instrument  phase  shifting  and  then  re- 
calculating the  cross-spectra. 
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After  all  the  required  power  and  cross-spectra  for 
a  block  haye  been  computed  and  stored,  the  sequence  is  repeated 
until  the  number  of  blocks  requested  haye  been  analyzed  or 
until  an  end-of-file  mark  is  sensed  on  the  coefficient  tape. 

Then  the  program  averages  the  spectral  estimates  both 
over  the  number  of  blocks  processed  and  over  the  analysis 
bandwidth  requested  in  the  input  data  cards.   The  standard 
deviation  from  each  mean  is  computed  in  a  similar  procedure 
from  the  formula 

E(Rn  -  Rp  1/2 


p(Rn  ■ 
=  l   N __ 

I   N-l 


where  N  is  the  number  of  samples  used.   Other  useful  information 
is  computed  and  printed  in  either  tabular  form  or  graphic  form 
as  specified  by  the  investigator. 

F.   IBM  360/67  TAPE  OPERATIONS 
1 .   Job  Control  Language 

Job  Control  Language  (JCL)  cards  provide  the  computer 
with  information  on  tape  mount  number,  tape  identification, 
disposition  of  tape;  identity  of  tape  data  file  to  process  and 
order  of  data  on  the  tape.   A  typical  JCL  card  used  for  tape 
processing  would  be: 

//GO.FT04F001  DD  UNIT=2400 ,VOL=SER=NPSl85 ,LABEL= ( 1 , SL) , 
DSNAME=MCKE01,DISP=( NEW, KEEP),  DCB  =  (DEN=2 
RECFM=VS,  BLKSIZE  =  8204) 
GO  -  This  group  indicates  data  is  going  to  be  classified 
under  the  subsequent  identify  parameters. 
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FT04  -  This  group  indicates  the  tape  is  to  be  placed  on 
logical  unit  (in  this  case  a  tape  mount)  four.   Also  the  two 
digit  number  (04)  must  match  the  number  specified  in  the  READ 
or  WRITE  instruction  in  the  program  (i.e.  READ  (4)  KMAX,  NCHAN, 
DATA).   The  JCL  does  not  determine  whether  "reading"  or 
"writing"  is  to  take  place.   The  program  statement  which  in- 
dicates the  logical  unit  (tape  mount)  also  indicates  the  pro- 
cess (READ  or  WRITE)  desired. 

F001  -  This  group  indicates  the  sequential  number  of 
passes  through  the  tape  which  has  occurred  up  to  this  point. 
In  this  case  it  is  the  first  pass. 

UNIT  =  2400  -  This  group  indicates  the  particular  unit 
will  be  a  nine-track  tape.   The  designation  for  a  seven-track 
tape  is  UNIT  =  2400  -1.  - 

VOL  =  SER  =  NPS  185  -  This  group  indicates  that  the  tape 
has  the  external  computer  center  marking  NPSI85.   If  a  seven- 
track  tape  were  used  it  could  have  a  user  specified  name 
(i.e.  MCKE01). 

LABEL  =  MCKE01  -  This  group  specified  the  Data  Set  Name 
(DSNAME)  of  the  data  specified  by  the  LABEL  group.   This  name 
may  be  changed  for  each  different  file  of  data;  however,  such 
a  procedure  is  quite  time  consuming.   It  is  much  quicker  to 
use  the  same  DSNAME  f or . a  whole  tape.   It  is  important  that 
when  ever  a  data  set  is  created  under  a  particular  name,  it 
must  always  be  referred  to  by 'the  same  name  when  reading 
from  the  tape. 
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DISP  =  (NEW,  KEEP)  -  This  group  specified  the  disposition 
of  the  tape  when  processing  (either  "reading"  or  "writing"). 
This  implies  it  is  a  newly  created  tape  (NEW)  and  it  is  to  be 
saved  (KEEP).   If  an  old  tape  were  to  be  read  and  the  data 
saved,  the  selection  would  be  DISP  =  (OLD,  KEEP). 

DCB  =  (DEN  =  2,  RECFM  =  VS ,  BLKSIZE  =  8204 ;  This  group  is 
the  Data  Control  Block.   The  most  often  used  density  (DEN)  is 
556  bytes  per  inch  which  is  denoted  by  2.   This  corresponds 
with  the  manual  setting  of  the  DENSITY  dial  on  the  face  of 
the  seven-track  tape  drive  unit .   The  tape  record  format 
(RECFM)  is  variable  spaced  (VS)  when  recording  and  digitizing 
on  the  seven-track  tape.   The  block  size  (BLKSIZE)  specifies 
the 'exact  number  of  bytes  in  one  block  of  data.   If  2048  words 
per  block  were  considered,  four  bytes  per  word  would  result 
in  8192  bytes  per  block.   If,  as  in  the  CONVERT  program,  two 
words  from  KMAX  and  one  word  from  NCHAN  give  an  additional  12 
bytes,  the  total  is  8204  bytes  per  block. 

2 .   Multi-file  Tape  Operations 

Most  frequently,  a  digital  tape  containing  several 
files  of  data  was  generated.   Under  normal  IBM  360  procedures, 
as  with  the  CONVERT  PROGRAM,  the  JCL  cards  are  used  to  indicate 
which  files  of  data  are  to  be  processed.   But,  in  the  FTOR  and 
SCOR  programs,  'JCL  cards  had  to  be  included  for  every  file, 
even  though  the  file  was  not  processed.   This  variation  to 
the  usual  JCL  procedures  was  necessary  because  of  the  SKIP- 
FILE  subroutine  in  both  programs  which  caused  the  files  that 
were  not  specified  for  analysis  to  be  skipped.  i 
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An  example  of  the  JCL  cards  for  the  normal  processing 
of  three  files  of  data;  where  it  was  desired  to  skip  the 
second  file  on  a  tape,  was  given  by  Jones  iRef.  5,  pg.  46]. 
3 .   Multi- Volume  Tape  Operations 

If  one  long  file  of  data  extends  onto  another  tape 
and  processing  of  this  long  data  file  is  desired,  multi-volume 
tape  procedure  becomes  quite  complicated  when  using  multi- 
tape programs,  its  use  should  be  avoided.   The  average  seven- 
and  nine-track  tapes  used  can  hold  more  than  1400  records  of 
2048  samples  per  record.   This  gives  a  capacity  for  more  than 
two  and  a  half  million  samples  per  tape. 

G.   PREPARATION  OF  CARDS  AND  TAPES  FOR  PSD  ANALYSIS 
•  1.   JCL  Cards  For  FTOR 

This  resume  of  JCL  techniques  for  using  the  FTOR 
program  are  taken  from  Jones  (1971).   The  following  procedure 
assumes  also  that  NPS  216  will  be  used  with  its  KMAX  set  to 
256.   If  changes  are  made  to  KMAX  and  the  program  is  taken 
from  another  tape  (as  discussed  under  Modifications  to  FTOR), 
the  same  argument  would  be  valid. 

The  FTOR  program  is  stored  on  NPS  216  in  File  1. 
The  program  is  also  storred  on  disk  (see  Wilson,  et .al . , 
(1969)).   To  use  FTOR  directly  from  the  tape,  the  following 
cards  should  be  used. 
JOB  CARD ~™r.„ __^.rw^^r^^T_^__ ,wr.w.... 

//FORT.  SYSIN  DD  UNIT  =  2400,' VOL  =  SER  -  NPS  216,  DISP  =  OLD, 
//       LABEL  =  (1,SL),DSN  =  UBCFTOR 
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JCL  CARDS  FOR  INPUT  &  OUTPUT  TAPES «-w 

//GO.  SYSIN  DD* 

—  CONTROL  CARDS  FOR  INPUT  &  OUTPUT  PARAMETERS ^— 

CBLANK  CARD) 
/*(See  Wilson,  et.  al.  iRef.  7,  Pg  33a  and  33bJ )  . 

If  the  FTOR  deck  is  used,  the  following  cards  should 
be  used: 

JOB * 

//  FORT.  SYSIN  DD* 

FTOR  PROGRAM  DECK 

JCL  CARDS  FOR  INPUT  &  OUTPUT  TAPES  

//GO.  SYSIN  DD* 

CONTROL  CARDS  FOR  INPUT  PARAMETERS  

(BLANK  CARD) 
/*  " 

Jones  [Ref.  5  p.  70]  gives  a  good  example  of  how  the 
input  cards  for  a  complete  FTOR  run  should  appear.   His  format 
assumes  that  the  FTOR  program  is  being  input  from  a  deck  of 
cards . 

,  The  calibration  card  referred  to  has  a  place  for  a 
primary  and  secondary  channel.   When  the  complete  signal  comes 
from  only  one  "primary"  input,  the  primary  channel  calibration 
card  is  the  only  one  used.   If  an  instrument,  such  as  the 
sonic  anemometer  is  used,  signals  need  to  be  added  yectorially 
to  give  the  actual  signal.   The  FTOR  program  allows  for  the 
input  of  two  separate  digital  signals.   The  calibration  cards 
permits  these  signals  to  be  added  vectorially  and  the 
spectrum  computed  for  a  single  signal.   The  alphameric 
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section  of  this  card  allows  for  signal,  parameter  identification, 
such  as  sampling  rate,  filter  setting,  analog  tape  channel 
number,  etc.  to  be  storred  with  the  Fourier  coefficients.   This 
alphameric  section  is  printed  as  the  title  on  all  output 
graphs . 

A  problem  might  arise  from  the  use  of  this  alphameric 
section  when  this  title  is  applied  to  long  data  lengths.   When 
short  sections  of  the  signal  are  analyzed  with  SCOR,  this 
title  information  is  printed  on  all  graphs.   Thus,  a  slight 
bookkeeping  problem  occurs  when  thirty  graphs  are  printed  with 
the  same  title.   Fortunately,  the  computer  center  delivers 
graphs  in  a  continuous  roll.   Identification  is  made  by  find- 
ing the  starting  point  and  then  numbering  them  sequentially 
from  that  point . 

As  mentioned  previously,  the  SKIPFILE  subroutine  in 
FTOR  made  it  necessary  to  include  JCL  cards  for  every  file 
on  the  nine-track  tape.   This  is  a  variation  to  the  normal 
IBM  360  procedure  in  which  the  JCL  card  specifies  which  file 
to  process.   The  SKIPFILE  routine  in  FTOR,  SCOR  and  FCPLOT 
allows'  for  the  internal  selection  of  data  files,  based  on  the 
parameters  on  the  data  cards  which  follow  the  tape  JCL  cards 
in  the  input  deck. 

2.   Modifications  of  FTOR  Program 

Jones  (1971)  found  it  was  necessary  to  make  several 
changes  so  that  FTOR  would  be  compatible  with  the  block  size 
specified  by  KMAX .   It  was  discovered  that  the  FTOR  program 
on  NPS  216  was  written  for  a  block  size  of  256  words.   The     ** 
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change  involved  setting  KMAX,  Card  18,  equal  to  the  number  of 
words  in  one  block  of  data.   Also  the  KMAX  dimension  statement, 
card  23,  was  similarly  changed. 

This  change  was  made  by  inserting  the  new  card  into 
the  program  deck.   Then,  in  order  that  the  whole  FTOR  deck 
need  not  be  input  each  time,  the  modified  program  was  written 
into  File  1  of  tape  NPS  223.   In  order  to  reduce  the  number 
of  different  tapes  used,  the  SCOR  program  was  written  into 
File  2  of  NPS  223.   If  too  many  different  tapes  are  used, 
confusion  arises  in  mounting  tapes.   Thus,  the  modified  FTOR 
and  SCOR  could  be  called  from  the  same  tape  by  the  use  of 
several  cards  rather  than  manipulating  an  entire  deck  of  cards. 
3.   JCL  For  SCOR 

The  cards  required  for  the  "running"  of  the  SCOR  pro- 
gram follow  the  same  order  as  those  required  for  FTOR.   The  - 
major  changes  are: 

LABEL  =  (2,SL)   The  program  is  stored  in  File  2 
DSN  =  UBCSCOR   The  data  set  name  is  UBCSCOR 
JCL  Cards  for  tapes  must  be  changed  to  read  the 
FTOR  coefficients  from  the  FTOR  generated  tape. 
Control  cards  for  analysis  must  be  appropriate  for 
SCOR  program,  (see  Wilson,  et_.al_.  (1969)  for 
details )  . 

Jones,  [Ref.  5  p.  71],  gives  a  good  example  of  how 
the  input  cards  for  a  complete  SCOR  run  should  appear.   His 
format  assumes  the  SCOR  program  is  input  from  NPS216.   This 
tape  can  be  used  without  the  KMAX  changes  needed  for  FTOR, 
because  the  input  to  the  SCOR  program  was  the  Fourier  co- 
efficient tape  which  had  a  uniform  format,  regardless  of  the 
block  size  used  in  the  nine-track  digital  input  tape. 
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The  selection  of  a  uniform  set  of  axes  should  be 
anticipated  so  that  graphic  outputs  can  be  overlaid  for  easy 
comparison  of  spectra.   If  logarithmic  bandwidths  (exponential) 
are  selected,  a  constant  bandwidth  bar  appears  on  the  log-log 
plot  due  to  the  exponential  nature  of  the  axes. 

Jones  iRef.  5 1 ,  points  out  that  the  SUBSEQUENT  ICMAX 
CARDS  referred  to  by  Wilson,  et_.  al.  [Ref.  7,  pg.  36]  in- 
dicate what  type  of  graphic  output  is  desired:  spectra,  cross- 
spectra  or  both.   The  spectrum  of  a  single  channel  is  specified 
as  if  it  is  a  cross-specturm  of  that  channel  with  itself.   For 
two  channels  (1  and  2),  a  spectrum  of  channel  one  and  a 
spectrum  of  channel  two  are  obtained  by  the  following  entries 
on  these  data  cards: 

.  First  card:    112 
Second  card:   2   2 

The  cross-spectrum  is  computed  only  is  each  channel  appears 

in  the  list  on  the  other  channels  card  as  seen  with  the  first 

card  above. 

4.   JCL  for  FCPLT 

This  program  is  used  to  get  a  Calcomp  plot  of  the 

amplitude  of  the  Fourier  coefficients  as  a  function  of  the 

frequency.   The  cards  required  for  "running"  the  program 

follow  the  same  order  as  required  by  FTOR  and  SCOR.   The  major 

changes  required  are:  . 

a.  LABEL  =  (3,SL)   This  program  is  stored  in  File  3- 

b.  DSN  =  UBCFCPLT  The  data  set  name  is  UBCFCPLT. 

c.  JCL  cards  for  tapes  may  be  the  same  as  those  used 
for  SCOR  if  the  same  data  files  are  to  be  used.  - 
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d.   Control  cards  for  analysis  parameters  required 
must  be  appropriate  Tor  FCPLT  as  specified  in  Wilson,  et  .al . 
iRef.  7,  pg-  59J. 

Jones  [Ref.  5,  pg.  72]  gives  an  excellent  example  of 
how  the  input  cards  for  a  complete  SCOR  run  should  appear. 
Figure  26  was  prepared  to  further  clarify  the  JCL  cards  needed 
for  each  program. 

5.   Suggestions  for  Efficient  Tape  Processing 

Since  the  techniques  involved  in  running  this  program 
can  be  quite  time-consuming,  much  attention  was  devoted  to 
the  problem  of  optimizing  the  number  of  runs  per  day.   Another 
problem  was  the  low  priority  for  these  programs.  (Priorities 
range  from  the  highest  priority,  class  A  to  the  lowest,  class 
K) .   The  low  priority  resulted  because  of  two  factors.   First, 
the  program  computation  time  (specified  on  the  green  JOB  card 
by  TIME  =  MM,SS,  where  M  is  minutes  and  S  is  seconds),  in  most 
cases,  took  more  than  four  minutes  of  Control  Processor  Unit 
(CPU)  time.   The  second  factor  was  that  the  program  required 
using  tape  input  and  output  units  which  are  automatically  put 
into  a  lower  class  because  of  the  time  required  to  mount  them 
and  access  the  information  on  them.   If  the  programs  and  the 
large  quantity  of  digitized  signals  could  be  completely 
stored  on  magnetic  disks,  higher  priorities  could  be  achieved. 
This  procedure  would  save  time  if  repeated  analyses  were  to 
be  conducted  on  a  data  set  which  did  not  change.   In  this 
study  many  different  signals  were  studied,  which  implied  the 
data  set  varied  from  one  sienal  to  another. 
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7-Track    Digitized    Data    Tape 


"/"/ CONVERT      PROGRAM      DECK 

//C0.FT02FO01    DD    UNIT=2400-1  ,UOL  =  SER=  MCKEfM  ,  LA-El  -  '  1 

//       DISP=(OLD,KEEP),DCB=(DEM=1,RECFI'=r,fcLK'-.IZE-r-ii'  ;' 

//GO.FT04F001  DD  UNI  T=  2400, V0L=SER=NPS1  P5  ,  LABEL-  (  1  '"■L1  ' 

DSNA«E«|V,CKE01,DISP=(NE!ll,KEEP),DCS=(DCN=2,«?EEFi»'=\.'S,eLKSIIEr  i\ 

Additional  JCL  cards  Included  for  eac*>  1 ...... 

File  for  which  conversion  is  desired. 


9-Track  Digitized  Data  Tape 


FTOR   PRUCRAM   DEC" 

0.FT)BF00O  DD  UN  I  T*  ?4nn,uTIL  =  SE  R  =  Nf  ,1  t'S  ,  DI  5P  =  0Ln  LAbEL-(l  Si) 
OSN*("CKEni,DCB=(OEN=?,RECFC=;S,SLKSIZE  =  ??n«)        ' 

Additional  JCL  cards  For  eac-  rilp  on  NPSiiS 

Because  of  SKIPFILE  routing  :l  card  'or  each 
file  on  mPSlPS  must  De  included. 
O.FT03F001   DD  UNI  T=2400,\/OL  =  SER=:\PS2??.  ,DI  3P=(tli£UI  KEFP)  LABEL-M  Sl^ 
DSN=C0Er-,DC8«(DEN=2,RECFM=y3,-LK5IZi>  :CA)      '  ' 

Additional  JCL  for  each  file  to  oe  written 

on  NPS223. 
O.SYSIN   DD  • 

Four  CONTROL  CARDS  for  each  file  ror  which 

Fourier  coefficients  are  to  be  computed. 
(  BLANK  CARD  TERMINATES  COMPUTATIONS  ) 


9-Track  Fourier  Coefficient  Tape 


//FORT.SVSIN  DD  UNI T= 2400, V0L=SERrNPS21 6 , DI SP=0LD. L AbEL* f 2. SL ) 
//       DSN=UBCSC0R 

//G0.FT03F001  DD  UNIT=2400,U0L  =  SER=NPS1B5, DISP= (OLD  ,  KEEP  )  LABEL=(1 
//      DSN=C0EF,DC8=(DEN=2,RECFMsyS,6LKSIZE=e204) 

Additional  JCL  cards  for  each  file  

on  NPS18S. Because  of  SKIPFILE  routine 


JCL  card  for  every 
must  be  Included. 

//GO.SYSIN  0D  • 

...POWER  SPECTRUM  PLOT  OF  TAPE  001 

-----------Three  cards  minimum  for 

which  spectrum  is  desired 
(  BLANK  CARD  TERMINATES 

/• 


file  on  NPS1P5 


j.d.mekendri; 

each  file  for- 


K.  .GCEAN0. , 


:0MPJTATI0NS  ) 


PSD  G^iplns       1>5D  Valuts 


Figure  26.   Program  Sequencing  and  JCL  Cards  Needed  for  PSD 
Analysis 
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a.   Program  Submission  Submission 

The  best  time  for  processing  digital  tapes  and 
debugging  particular  problems  was  found  to  be  during  the  summer 
and  winter  break  when  the  work  load  at  the  computing  center 
was  very  light.   During  that  time,  generally,  about  ten  runs 
per  day,  depending  on  the  length  of  the  JOB,  could  be  achieved. 
If  all  was  going  well  only  one  run  per  day  was  often  adequate; 
however,  when  problems  developed,  an  immediate  debugging  run 
was  essential.   It  was  found  that  oversights  on  the  part  of 
the  programmer,  in  the  form  of  errors  in  JCL  and  control  cards, 
the  "dumb"  computer  was  only  doing  what  it  had  been  instructed 
to  do.   However,  occassional  computer  malfunctions  did  occur. 

The  next  best  time  for  program  processing  was  on 
weekends  at  the  beginning  of  the  quarter.   As  the  quarter  pro- 
gressed, the  weekend  would  offer  several  runs  per  day.  Toward 
the  end  of  the  quarter,  due  to  heavy  work  load,  only  one  run1 
at  around  3  A.M.  resulted,  if  the  job  had  been  input  prior  to 
10  A.M.  of  the  previous  day.   In  general,  turn  around  averaged 
24  hours. 

The  most  CPU  used  for  any  program  sequence  was 
ten  minutes.   Due  to  the  fact  that  each  analysis  uses  different 
parameters,  the  only  rule  of  thumb  which  was  found  useful  was 
that  the  FTOR  procedure  took  about  22.5  seconds  per  record  (204  8 
samples  per  record)  to  compute  the  Fourier  coefficients.  Thus, 
400  records  required  490  seconds  of  CPU  time.   The  FTOR  pro- 
gram was  the  most  time  consuming  thus  this  figure  (22.5  seconds) 
can  be  used  as  an  upper  limit  in  estimating  CPU  time  for 
other  programs . 

76 


b.   Stacking  Programs 

Another  technique  used  to  affect  faster  PSD 
analysis  involved  submitting  CONVERT,  ETOR  and  SCOR  under  one 
JOB  card.   The  reason  the  programs  were  "stacked"  was  to  per- 
mit CONVERT,  FTOR  and  SCOR  programs  to  run  sequentially. 
Btherwise,  it  was  necessary  to  get  the  results  of  each  program 
before  the  next  program  could  be  input .   If  no  mistakes  were 
made  in  the  JCL  cards  and  if  no  problems  existed  with  the 
tapes,  the  complete  analysis  would  be  made  in  one  run.   If  a 
problem  in  any  step  was  encountered,  corrective  action  was 
taken  and  the  remaining  programs  were  run  individually.   The 
option  as  to  when  to  stack  the  programs  and  when  to  run  them 
individually  varied  with  the  number  of  files  on  the  tapes. 
The  higher  the  number  of  files,  the  more  JCL  and  control 
cards  required,  and  the  greater  the  chance  for  errors. 
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IV.   EXPERIMENTAL  PROCEDURE 

The  initial  plan  for  the  investigation  of  the  noise  pro- 
blem,  involved  digitization  and  spectral-analysis  of  real 
signals.   These  were  to  be  sine  waves,  square  waves,  ramps, 
and  random  noise.   Since  these  input-signal  characteristics 
were  readily  known,  the  Fourier-transform  and,  the  power 
sectra  were  known  functions.   The  first  signal,  a  10HZ  sine 
wave,  did  not  give  the  expected  spectral  values.   The  digital 
program  was  checked  to  see  if  its  calculations  were  valid  and 
then  the  A/D  procedure  was  checked.   For  purposes  of  clarity, 
a  chronological  discussion  of  the  experimental  procedure  will 
follow. 

A.   ANALYSIS  OF  PURE  SIGNALS 

Though  the  first  series  of  10HZ  sine  waves  showed  an 
expected  energy  peak  at  10HZ,  the  exponential  decrease  of 
energy  with  increasing  frequency  was  not  expected  for  a  sine 
function.   An  assessment  was  made  that  the  problem  could  be 
either  in  the  actual  A/D  step  or  in  the  spectral-analysis 
programs.   Before  definite  conclusions  concerning  this  pro- 
blem and  the  noise  sources  could  be  made,  it  was  necessary 
to  gather  baseline  information  on  computer  analysis  of 
theoretically  pure  signals. 

1 .   Computer  Generated  Digital  Sine  Function 

The  program  listed  in  Appendix  I  was  used  to  generate 
a  simulated  digital  sine  signal  and  the  computed  samples 
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were  stored  on  a  digital  nine-track  tape.   The  format  of 
the  tape  was  made  to  be  compatable  for  input  of  this  data  to 
the  FTOR  program.   In  the  sine  generation  program,  the  peak 
voltage  could  be  varied  by  changing  the  sine  amplitude,  and 
the  the  sampling  rate  could  be  changed  by  varying  At.   A 
standard  block  size  of  2048  samples  per  block  was  maintained 
throughout  this  study.   The  CONVERT  procedure  was  not  needed 
because  the  data  was  already  in  a  hexidecimal  format. 
2 .   Computer-Generated  Digital,  Random  Signal 

The  next  step  was  to  test  the  FTOR  and  SCOR  spectral 
density  analysis  of  a  signal  with  a  wide  range  of  frequencies. 
To  do  this,  a  Gaussian  signal,  which  has  a  flat  spectral 
density  function,  was  used.   The  program  listed  in  Appendix 
II  generated  a  simulated  digital  random  signal, by  "calling" 
the  computer  sub-routing  RANDU .   The  peak  voltage  values 
could  be  changed  by  altering  the  constant  multiplier  of  YFL; 
the  sampling  rate  could  be  changed  by  altering  the  sampling 
rate  specified  on  the  FTOR  input  data  deck.   The  peak  amplitude 
was  maintained  at  10  volts;  however,  the  different  sampling 
rates  were  investigated.   Since  the  data  was  stored  on  a  nine- 
track  tape  with  a  format  compatible  with  FTOR,  the  CONVERT 
program  was  not  used. 

B.   A/D  CONVERSION  AND  PSD  OF  LABORATORY  SIGNALS 

Once  characteristics  had  been  established  for  computer 
processing  of  pure  control  signals,  the  next  step  was  to 
digitize  actual  signals.   It  was  decided  that  a  random  or 
Gaussian  signal  would  give  all  the  information  required,  and 
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thus,  the  sine,  ramp,  and  square  wave  signals  could  be  by- 
passed.  The  random  signal  was  expected  to  give  optimum  power- 
spectral  density  information  for  purposes  of  the  study. 
1 .   Random  Signal 

a.  Single  Channel  Digitization 

An  Elgenco  model  603  A,  Gaussian  noise  generator 
was  used  to  give  a  random  noise  output.   Its  characteristics 
are  listed  in  Appendix  III.   A  frequency  setting  of  5Hz  to 
20  KHz,  attenuation  schale  XI . 0 ,  output  voltage  reading  2.62 
Vrms  was  input  to  a  Khron-Hite  filter  model  3321  set  at  2000 
Hz,  low  pass  max.  flat  and  Odb  gain.   The  filter  output  was 
input  to  the  Hybrid  computer  through  an  operational  amplifier 
with  a  ten  volt  gain.   A  sampling  rate  of  5000  SPS  was  selected, 
A  total  of  4l  seconds  of  the  signal  was  digitized  onto  a 
seven-track  tape.   An  End-of-File  mark  was  written  onto  the 
tape  by  typing  the  EOF  option  on  the  teletype  keyboard. 

b.  Dual  Channel  Digitization 

The  second  file  on  the  tape  mentioned  above  was 
filled  with  data  from  two  input  channels  which  were  digitized 
simultaneously.   The  Elgenco  noise  generator  described  above 
was  used  as  one  input,  and  the  other  input  was  from  the  random- 
noise  generator  which  is  built  into  the  Ci  5000. 

The  noise-generator  was  utilized  with  the  same 
settings  as  above;  however,  the  filter  was  reset  to  1880Hz. 
The  random-noise  generator  from  the  Ci  500Q  was  input  to  a 
Krohn-Hite  filter  also  set  at  1880  Hz.   The  output  of  this 
filter  was  fed  into  a  different  operational  amplifier  with 
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a  gain  of  10.   The  sampling  rate  was  lowered  to  4000  SpS. 

A  total  of  25.5  seconds  of  signal  was  digitized.   Though  the 

single  and  dual  channel  cases  both  considered  100  records, 

this  shorter  digitizing  time  occurred  because,  two  samples 

were  being  taken,  simultaneously,  at  the  rate  of  4000  SPS: 

2048  Samp  x  100  Blks 

4000  Samp      ~      =   25*5  Sec * 
Sec    x      d 

An  End-of-File  mark  was  again  written  on  the  tape  to  end  this 

file  of  data. 

The  CONVERT  program  as  used  to  convert  from  seven- 

to  nine-track  tape.   The  program  SCOR  allowed  the  computation 

and  plotting  of  the  spectrum  of  each  channel  and  the  eo- 

spectrum  and  the  quad-spectrum  of  one  channel  with  another. 

2 .   Random  Signal  and  Sine  Signal 

To  determine  whether  a  sine  wave  could  be  picked  out 
of  the  random  noise,  a  1000  Hz  sine  and  random  signal  were 
digitized.   An  attempt  to  digitize  a  single  channel  of  the 
sine  combined  with  the  random  signal  failed  due  to  a  faulty 
patch  on  the  analog  board.   PSD  values  were  obtained,  seemingly 
because  the  open  amplifier  actually  picked  up  stray  signal.  The 
sine  amplitude  was  increased  and  the  a  second  file  was 
digitized.   The  signals  in  these  two  files  were  digitized  at 
5000  SPS  and  filtered  at  2000  Hz. 

Dual  channel  digitized  samples  of  the  1000  Hz  sine 
with  amplitude   ±20  volts  and  Gaussian  signals  filled  the 
third  file.   The  amplitued  was  increased  to  ±30  volts  and 
the  two  separate  sine  and  Gassian  signals  digitized  into  the 
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fourth  file.   The  signals  in  these  third  and  fourth  files 
were  digitized  at  4000  SPS  and  filtered  at  2000  HZ. 

C.   DATA  FROM  GEOPHYSICAL  SIGNALS 

Atoraspheric-temerature  and  velocity  signals  have  pre- 
viously been  recorded  on  one-inch  magnetic  tape  by  Boston 
[Ref.  1].  These  signals  were  used  as  a  final  check  on  the 
system  to  determine  if  correct  spectral  values  could  be 
achieved. 

The  signals  were  reproduced  on  a  Sangano  Model  3562  FM 
tape  recorder  at  60  ips.   The  tape  playback  output  was  filtered 
at  1000  H'z  for  the  temperature  signal,  and  at  2000  Hz  for 
velocity,  the  differentiated  velocity  and  temperature  signals. 
The  sampling  rate  for  the  temperature-signal  digitization 
was  2000  SPS  and  the  other  three  signals  were  sampled  at  4000 
SPS.   The  filter  setting  was  low-pass  max,  flat,  Odb  gain. 


82 


V.   ANALYSIS  OF  RESULTS 

A.   PSD  OF  COMPUTER  GENERATED  SIGNALS 

1 .  Sine  Wave 

Figure  27  was  the  spectral  plot  of  a  computer- 
generated  sine  wave.  The  PSD  values  were  computed  for  24  re- 
cords of  signal  giving  a  total  signal  length  of  11.9  seconds. 
The  total  integrated  power  was  .499  V  ;  However,  the  power  in 
the  8.12  Hz  band  centered  about  7.11  Hz  had  a  total  of  .495  V2 
(8.13  Hz  X  6.09  X  10"2V2/Hz).   Essentially,  all  the  power  was 
contained  in  the  band  between  2.05  Hz  and  11.17  Hz.   Though 
this  appeared  to  be  a  wide  band,  the  whole  region  from  11.7  Hz 
to  263.1  Hz  had  less  than  .8  percent  of  the  total  power: 

•4".499495  x  100  -,8* 

Another  test  run  with  a  200  Hz  sine  wave  with  a 
peak-to-peak  amplitude  of  10  volts,  proved  inconclusive  due 
to  an  error  in  specifying  the  sampling  rate  on  the  FTOR  data- 
card input.    A  sampling  rate  of  400  SPS  was  specified,  rather 
than  500  SPS,  which  was  the  actual  rate  used.   The  expected 
total  power  value  of  50. V"-  was  achieved;  however,  the  error 
with  the  sampling  rate  shifted  the  frequency  peak  from  200  Hz 
to  159  Hz.   Though  the  error  produced  erroneous  results,  the 
effect  of  not  specifying  the  correct  sampling  rate  was  observed 

2 .  Gaussian  Noise 

Figure  28  was  the  PSD  plot  of  the  computer-generated 
random  signal.   The  digital  samples  of  the  signal  simulated 
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a  random  signal  of  9.8  seconds  in  length,  sampled  at  5000  SPS. 

The  spectral  level  was  found  to  be  very  flat  with  increasing 

-3   2 
frequency.   Average  spectral  density  of  3.30  X   10   V   /Hz 

-3   2 
varied  from  a  high  value  of  3 .45  X  10   V  /Hz  to  a  low  value 

of  3.25  X  10""-*  V  /Hz .   Actual  variance  was  very  low  and  quite 

uniform.   No  frequency  spikes  were  observed  and  the  conclusion 

was  drawn  that  the  random  number  generating  sub-routine  RANDU 

produced  a  true  random  series  for  at  least  the  first  50,000 

numbers.   It  was  noted  that  the  random-number  generating 

39 
capacity  of  this  sub-routine  is  2    numbers  before  the  series 

17 
repeats  itself.   Since  only  2  '  numbers  were  used,  the  full 

potential  of  the  number  generator  was  not  fully  tested. 

Figure  29  is  the  SPD  plot  of  the  same  random  series 

sampled  at  1000  SPS  rather  than  5000  SPS.   This  changed  the 

record  length  to  2.05  seconds,  and  since  the  results  for 

100  records  was  computed,  the  total  length  of  signal  sampled 

was  205  seconds.  Although  the  sampling  rate-change  produced  a 

higher  PSD  value,  the  mean  showed  little  variation.   The  mean 

"2   2  ? 

level  was  about  1.57  X  10   V  /Hz  with  a  low  of  1.62  X  10"^ 

2  2 

V  /Hz  and  a  high  of  1.72  V  /Hz.   The  variance  around  the  mean 

level  was  very  small.   As  expected,  averaging  over  a  fewer 
number  of  records  (24  versus  100),  the  variance  was  higher. 
B.   PSD  LABORATORY  SIGNALS 
1 .   Single  Channel  Sine 

a.   Signal  Leakage  Into  Open  Amplifier 

Figure  30  was  the  PSD  plot  of  the  Spectrum  of 
signals  with  peak  voltages  of  ±2.0  and  ±3.0.   The  signals 
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Figure  30.   Signal  Leakage  into  Open  Amplifier 
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were  picked  up  by  an  open  input  amplifier,  which  was  next  to 

the  amplifier  into  which  the  signals  were  actually  input.   The 

open  amplifier,  whose  output  was  being  digitized,  acted  like 

an  antenna  in  picking  up  these  stray  signals.   The  plots  show 

considerable  consistency  and  several  conclusions  can  be  made. 

The  spectral  peak  of  the  ±3  volt  signal  was  higher 

than  that  of  the  ±2  volt  signals.   The  peak  PSD  was  2.68  X  10~5 

2 
V  /Hz  for  the  lower  and  6.58  X  10"5  V2/Hz  for  the  higher  signal 

These  values  multiplied  by  the  band  width,  78. 1  Hz  in  both 
cases,  gave  power  of  2.09  X  10"3  V2  and  5.12  X  10~3  V2 
respectively.   The  power  ratio  produced  by  this  difference  in 
output  voltage  was 

Power  Ratio  =  6.56  X  10~5 

2.68  X  10-3  -  2.45. 

Since  this  was  a  power  ratio,  the  voltage  ratio"  is  the  square 
root  of  the  power  ration,  or 

Voltage  Ratio  =   2.45  =  I.56 

The  actual  voltage  input  into  the  open  amplifier 

was  computed  from  the  power  of  each  signal.   The  lower  power 

-3  2 
2.09  X  10   V  resulted  from  an  input  voltage  of  ±4.5  X  10~2V, 

~3  2 
and  5.12  X  10  V  resulted  from  an  input  of  ±7.2  X  10~2V. 

The  expected  voltage  ratio  was  computed  using  the 
observed  input  voltage  values: 

Voltage  Ratio  =  3.0 

2.0    1,& 

The  power  ratio  was; 

Power  Ratio  =  (1.5)2  =  2.25 
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The  ratio  of  the  signal  voltage  in  the  input 
line  to  the  actual  computed  signal  voltage  gave  the  actual 

percent  of  signal  picked  up  by  the  open  amplifier. 

_2 
4^2^Q-1Q —  X   100  =  2.3$  for  the   2V  signal 



7,2 X   10 —  X  100  =  2.4$  for  the   3V  signal 

Assuming  the  leakage  was  coming  from  the  voltages 
in  the  input  lead  to  the  open  amplifier,  the  percent  of 
signal  leakage  would  be  the  ratio  of  the  voltage  in  the  input 
lead  to  the  signal  voltage  computed  from  the  PSD.   For  the 
±2  V  and  ±3  V  signals  the  percent  of  signal  picked  up  by  the 
open  amplifier  was  2.3  percent  and  2.4  percent  respectively: 

4.6  X  10~2 


2.0 


X  100  =  2.3$ 


7t2  X  10~    X  100  =  2A% 
j  .  u 

If  the  signal  was  leaking  from  the  closed  amplifier, 

the  leakage  was  .2  percent  for  both  cases. 

4.6  X  10~2 


2.0 


x  ioo  =  .2: 


_2 

7-2  x  1Q     X  100  =  .2% 
3.0 

Thus,  only  a  small  signal  leakage  was  observed  and 
it  was  independent  of  signal  amplitude. 

b.   Effect  of  Increasing  Signal  Amplitude 

Figure  31  is  a  PSD  plot  showing  the  effect  of  in- 
creasing the  voltage  of  the  signal  going  into  the  data  taking 
amplifier.   The  input  signals  had  peak-to-peak  voltages  of 
±20  volts  and  ±30  volts. 
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The  PSD  for  the  peaks  3-92  V2/Hz  and  9.29  V2/Hz 

were  computed  for  a  bandwith  of  68.4  Hz.   The  power  was 

268V  and  635V  respectively.   This  gave  a  power  ratio  of; 

Power  Ratio  =9.29 

3732  =  2.37 

The  voltage  ratio  is  the  square  root  of  the  power  ratio  or 
Voltage  Ratio  =  *\  2.37  =  1.54 

Using  the  power  spectral  density  to  compute  the 

o  2 

voltage  ratio,  268V^  implied  an  input  of  16.4  V    and  635  V 

rms 

implied  an  input  of  25-2  V    .   This  implied  peak-to  peak 

rms 

voltages  of  ±23. 2V  and  ±35. 6V.   Due  to  the  inaccuracy  involved 

in  reading  peak-to-peak  voltages  from  the  oscilloscope  on  the 

Ci  5000,  the  observed  inputs  of  ±20V  and  ±30V  could  have  been 

±5V  in  error. 

Assuming,  the  inputs  were  of  ±20V  and  ±30V,  the 

expected  voltage  ratio  would  have  been: 

Voltage  Ratio  =  30  -  ]  5 

20 

and  the  expected  power  ratio  would  have  been: 

Power  Ration  =  (1.50)2  =  2.25 

The  observed  and  computed  power  ratios  compared 
favorable,  and  the  difference  between  ovserved  and  computed 
peak-to-peak  voltages  (  20V  Vs.  ±23. 2V  and  ±35. 6V)  were  with- 
in acceptable  limits. 

The  theoretical  power  for  sine  waves  of  ±20V  and 
±30V  was  found  from  the  formula; 
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p 

where  V  is  peak-to-peak  voltage.   This  gave  power  of  200  V 

2 

and  450V  for  the  two  voltage  signals  respectively.   The 

difference  between  expected  power  level  and  the  power  level 
derived  from  the  spectral  plots  was  assumed  to  be  due  to  the 
error  in  reading  the  input  signal  amplitudes. 
2 .   Single  Channel  Gaussian  Signal 

Figure  32  was  the  spectral  plot  of  40. 06  seconds  of 
a  random  signal  sampled  at  5,000  SPS .   The  spectrum  level 
was  very  flat  to  about  1.5  KHZ.   Beyond  1.5  KHz,  rapid  decrease 
in  the  power  spectral  density  with  increasing  frequency  was 
to  be  expected.   The  3db  down  point  occurred  at  2.0  KHZ.   The 
slope  of  the  filter,  as  specified  in  the  equipment  characteristics, 
was  -48db/octave .   The  observed  slope  was  very  close  to 
-96db/octave .   This  value  would  be  expected  if  two  filters  had 
been  cascaded,  but  this  was  not  the  case.   The  spectrum  was 
quite  free  of  noise  spikes  and  had  no  60  Hz  harmonics  present. 

If  60  Hz  noise  was  present,  its  level  was  well  below  -17.5db/Hz. 

-2  ? 
The  spectral  level  of  1.73  X  10  V  /Hz  compared 

favorably  with  the  input  signal.   Specifications  for  the 

Elgenco  noise  generator  gave  a  spectral  density  of  approximately 

-3 

5  X  10   V/   Hz  at  IV    .   The  input  signal  had  a  meter  read- 

rms 

ing  of  2.62V    .   The  gain  factor  was  10.   The  computed  spectral 
&         rms       & 

density  was  2.69  X  10~2V2/Hz. 

/5  x   1Q-3V\2  00  _?   ? 

\      Hz  /  (2.6r(10r    =    2.69   x   10    "TVHz 

Since  no  accurate  spectral  density  information  on  the  noise 
generation  was  available,  these  values  are  considered  to 
compare  favorably  . 
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3 .   Two  Channel  PSD  of  Gaussian  Signals 

Figure  33  is  the  PSD  plot  of  one  of  two  channels  of 
random  noise  which  were  digitized  simultaneously.   One  spectrum 
resulted  from  the  Elgenco  Gaussian  noise  generator  and  the 
other  spectrum  resulted  from  the  high  frequency  random-noise 
generator  of  the  Hybrid  computer.   The  flat  spectrum  of  the 
Elgenco  generator  was  contrasted  by  the  very  "colored"  spectrum 
of  the  Ci  5000  random-noise  generator,  whose  low  frequencies 
contain  much  power.   Because  of  its  deviation  from  the  Gaussian 
characteristics,  future  use  of  the  Hybrid  noise-generator 
should  be  avoided  when  a  Gaussian  generator  is  wanted. 

C.   PSD  OF  TURBULENCE  SIGNALS 

1 .   General  Signal  Characteristics  Found;  Comparison 

with  Previous  Results 

a.   Temperature  Signal 

Figure  3^  showed  the  PSD  of  the  temperature  signal 
recorded  by  Boston  [Ref.l  ].   A  comparison  was  made  between  the 
PSD  valued  obtained  from  56  seconds  of  signal  that  was 
digitized  and  analysed  at  the  Naval  Postgraduate  School  Facility, 
and  PSD  the  values  obtained  by  Boston  from  the  same  section 
of  signal  that  was  digitized  and  analyzed  at  the  University  of 
British  Columbia  facility. 

The  general  PSD  characteristics  from  both  analysis 
compared  favorably  in  slope  magnitude  and  60  Hz  harmonic  peaks. 
The  strong  -5/3  slope  region  was  evident  on  both  spectra  be- 
tween 10  Hz  and  120  Hz.   Significant  harmonic  levels  were 
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evident  in  the  Naval  Postgraduate  School  study  at  the  first 
(60Hz)  and  third  (l80Hz)  harmonics  and  a  less  significant 
ninth   (5^0  Hz)  harmonic.   Boston's  results  showed  a  significant 
first  harmonic.   The  third  harmonic,  though  detectable,  showed 
a  very  minor  level;  the  ninth  harmonic  was  not  detected.   A 
slight  increase  in  the  slope  to  approximately  -7/3  was  ob- 
served between  200  Hz  and  400  Hz  (600  Hz  in  Boston's  results) 
followed  by  a  marked  decrease  in  the  slope  between  600  Hz  and 
1000  Hz. 

The  major  difference  between  the  two  spectra  is 
the  constant  power  level  difference.   The  Naval  Postgraduate 
School  results  were  higher  by  a  constant  factor  of  400  which 
implied  that  a  voltage  gain  difference  of  20  was  present  in 
the  Naval  Postgraduate  analysis. 

As  can  be  seen  from  Figure  35,  in  which  each 
Naval  Postgraduate  School  PSD  value  has  been  reduced  by.--a 
factor  of  400,  the  spectra  for  the  two  analyses  compare  vary 
favorably  in  slope  and  relative  magnitude.  A  temporal  PSD 
plot  described  in  a  later  section,  would  have  been  very  help- 
ful in  detecting  noise  in  Boston's  results;  however,  time  did 
not  permit  its  production. 

b.   Differentiated  Temperature  Signal 

A  comparison  of  short  duplicate  sections  of  the 
differentiated  temperature  signal  was   not  undertaken  in  this 
study.   A  long  section  of  the  signal  was,  however,  analysed. 
Figure  36  shows  the  comparison  for  the  PSD  from  five  minutes 
of  signal  obtained  in  this  study  with  the  PSD  from  three 
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different  sections  of  data  obtained  by  Boston  (Ref.'lJ.   The 
general  characterisitics  were  observed  to  agree  favorable, 
though  the  Naval  Postgraduate  School  results  were  greater  by 
a  constant  factor  of  300..  In  Figure  37  the  Naval  Postgraduate 
School  results  were  reduced  by  a  factor  of  300  and  were  found 
to  follow  closely  the  results  from  E  l'(B)  and  E  l'(C).   The 
first  harmonic  of  60  Hz  was  found  in  the  Naval  Postgraduate 
School  and  University  of  British  Columbia  data;  however,  the 
strong  third  harmonic  found  in  the  NPS  results  .wasn 't  evident 
in  the  UBC  results. 

c.   Velocity  Signal 

Figure  38  shows  the  PSD  of  a  section  from  the 
velocity  signal  recorded  by  Boston  on  tape  203  (1).  It  com- 
pares with  a  section  analysed  by  Boston  and  is  referred  to  as 
U(A).   The  original  analysis  was  conducted  for  a  sample  which 
was  about  56  seconds  in  length.  The  original  sampling  rate  was 
200  SPS  which  gave  a  Nyquist  frequency  of  l.OKHz.  The  Naval 
Postgraduate  School  analysis  was  conducted  on  the  same  section 
of  signal,  using  a  sampling  rate  of  4000  SPS  and  a  block  size 
of  2048  samples  per  block.   Since  only  56  records  were  analyzed, 
the  total  length  of  signal  was  only  51-2  seconds.   Though  the 
Naval  Postgraduate  School  results  are  based  on  a  slightly 
shorter  signal,  the  PSD  plots  show  quite  similar  results. 
The  power  level  of  the  University  of  British 
Columbia  analysis  is  lower  than  the  Naval  Postgraduate  School 
analysis  by  a  factor  of  340,  which  implies  a  voltage  gain 
difference  of  18.4  existed  between  the  two  sets  of  results. 
Harmonics  of  60  Hz  were  found  in  both  cases  in  that  the  first 
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and  third  harmonic  peaks  stood  well  above  the  -5/3  slope  of 
the  signal.   A  marked  deviation  in  results  occurred  in  the 
PSD  plots  for  frequencies  higher  than  250  Hz.   The  University 
of  British  Columbia  results  showed  a  spectral  peak  at  about 
320  Hz  whereas  the  Naval  Postgraduate  School  results  showed 
a  peak  at  620  Hz.   Since  the  Nyquist  frequency  of  this  study 
was  1000  Hz  higher  than  Boston's  results  for  U(A),  previous 
values  had  not  been  reported  for  the  frequency  range  1000  to 
2000  Hz.   Indications  from  Figure  39  are  that  the  slope  has 
decreased  from  the  -5/3  value  to  a  value  of  -2.5/3. 
2 .   New  Results  Obtained 

Based  on  the  encouraging  comparison  of  results  obtained 
in  this  study  with  the  results  obtained  by  Boston  at  University 
of  British  Columbia,  the  values  obtained  for  the  longer  time 
period  of  five  minutes  appear  to  offer  new  insight  into  the 
statistical  properties  of  the  geophysical  processes  measured. 
Generally,  the  results  obtained  indicate  that  the  statistical 
properties,  of  judiciously  chosen  short  sections  of  a  signal, 
do  give  valid  indications  of  the  nature  of  over-all  processes 
under  consideration- 

With  the  increased  capability  afforded  by  computer 
analysis  of  data,  new  insights  may  be  achieved  toward  viewing 
natural  phenomena. 

a.   Temporal  Variations  in  the  PSD 

A  different,  but  not  by  no  means  new  way  of  dis- 
playing computer  PSD  values  of  geophysical  processes,  is 
through  temporal  PSD  analysis.   This  technique  allows  the 
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investigator  the  opportunity  to  relate  fluctuations  of  the  PSD 
with  signal  fluctuations.   A  section  of  a  signal  which  appeared 
to  have  a  nominal  signal  amplitude,  as  in  Figure  4la  had  a 
low  PSD  value  (as  seen  in  the  section  of  PSD  plot  marked  10.24 
in  Figure  42.   Though  the  characteristics  of  the  PSD  are  in- 
deed determined  by  sampling  rate,  the  number  of  samples  in  each 
digitized  record,  filter  settings,  background  noise,  etc., 
meaningful  results  can  be  achieved  within  these  limitations. 

Since  the  identity  of  each  digitized  record 
(digitized  block  of  data)  was  maintained  throughout  the  analysis 
under  UBC  FTOR  and  UBC  SCOR,  the  PSD  of  sequential  sections 
of  the  signal  could  be  computed.   Temporal  variations  in  the 
PSD  of  temperature  and  velocity  signals  were  drawn  from  the 
results  obtained  by  sequentially  analyzing  a  constant  number 
of  records.   Plots  were  also  drawn  from  the  results  of 
analyzing  an  increasing  number  of  records,  beginning  with  the 
first  record.   These  plots  showed  that  the  fluctuating  PSD 
became  more  stationary  when  more  and  more  data  is  analyzed. 
(1)   Temperature. 

Figure  42  shows  the  temporal  variation  of 
the  temperature  PSD.   A  total  of  300  records  were  analyzed  to 
give  a  total  time  of  307.2  seconds.   To  compute  the  PSD  values, 
thirty  passes  were  made  through  the  Fourier  coefficients  with 
the  PSD  being  computed  for  each  set  of  ten  records.   The  time 
difference  between. each  PSD  was  10.24  seconds. 

Figure  43  showed  how  the  temporal  variations 
in  the  PSD  are  smoothed  out  by  taking  successively  longer 
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records  of  data.   The  first  PSD  was  computed  from  10  records, 
starting  with  record  number  one.   The  second  PSD  along  the 
time  axis  was  computed  from  20  records,  starting  with  the 
first  and  so  forth.   The  last  PSD  gave  the  average  of  a  total 
of  102.4  seconds  of  temperature  fluctuations,  which  was  com- 
puted from  100  records.   The  smoothing  effect  was  quite  evident 
from  the  fact  that  the  first  PSD  was  a  much  lower  value, 
though,  its  effect  was  quickly  smoothed  by  the  averaging  pro- 
cess . 

(2)   Velocity. 

Figure  44  shows  the  temporal  variation  of 
the  velocity  PSD.   A  total  of  300  records  were  analyzed.  Be- 
cause of  the  higher  sampling  rate  of  4000  SPS,  the  total 
length  of  signal  shown  is  only  153-7  seconds: 

SSjgj  Samp/Record    30o  Records 
4000  Samp/Sec 

Thirty  separate  PSD  values  were  computed  by  sequentially 

analyzing  10  records  at  a  time.   The  time  difference  between 

each  PSD  curve  was  5.12  seconds. 

Figure  45  showed  the  smoothing  effect 

achieved  by  analysing  an  increasing  number  of  records  of  the 

velocity  signal.   The  overall  trend  in  smoothing  was  not  as 

apparent  for  this  signal  due  to  the  fact  that  the  original 

section  of  signal  analyzed  gave  a  good  statistical  description 

of  the  process. 

b.   PSD  Five  Minutes  of  Signal 

Figure  46  showed  the  temperature  PSD  values  for  a 

long  length  of  data  (5  minutes)  compared  with  data  from  a 
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short  section  of  signal.   The  two  lengths  compare  very  closely, 
implying  that  for  this  situation,  the  shorter  length  of  a 
minute  would  haye  giyen  statistics  representative  of  longer 
sections,  (for  the  frequency  range  1  Hz  to  lKHz). 

Figure  47  showed  the  slope  characteristics  for 
the  slong  record  length  (5  minutes).   A  definite  increase 
(-7/3)  in  the  slope  was  noted  from  about  70  Hz  to  about  300  Hz. 
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VI.   CONCLUSIONS  AND  RECOMMENDATIONS 

A.  CONCLUSIONS 

The  following  conclusions  were  reached  concerning  the 
problem  of  the  possible  input  of  noise  into  the  Analog-to- 
Digital  conversion  procedure  and  into  the  PSD  analysis  pro- 
cedure with  signal  inputs  of  real  data. 

1 .  PSD  Programs 

No  significant  noise  sources  were  found  to  exist  with- 
in the  computational  program  of  the  Naval  Postgraduate  School 
FFT  package.   The  PSD  computations  for  noise  -free  sine  and 
random  signals  agreed  very  well  with  theory. 

The  programs  FTOR,  SCOR  and  FCPLT  were  found  to  be 
excellent  for  obtaining  power  spectra  from  large  quantities 
of  time  series  data. 

In  the  CONVERT  program  used  in  a  previous  study,  it 
was  found  that  a  factor  was  missing  which  changed  the  octal 
base  number  to  the  hexidecimal.   A  corrected  version  of  the 
program  was  used  and  PSD  results  indicate  the  procedure  is 
functioning  correctly. 

PSD  values  from  four  atmospheric  turbulence  signals 
which  were  compared  with  results  obtained  from  other 
computational  facilities  showed  very  close  correlation. 

2.  Analog-to-Digital_  Conversion 

No  significant  noise  sources  were  found  to  be  inter- 
fering with  the  digitization  procedure  carried  out  on  the 
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Hybrid  Computer.  Patchboard  noise  which  had  previously  been 
reported  as  excessive  by  Jones ,  £Ref.  5J  was  noted;  however, 
its  presence  in  the  PSD  plots  of  sine  and  random  signals 

were  not  detected.   Even  the  noise  picked  up  by  open  10  gain 

-6   p   -1 
amplifier  had  a  very  low  level,  on  the  order  of  10   V   Hz 

Early  in  this  study  it  was  discovered  that  the 
digitization  procedure  was  missing  several  data  samples  at  the 
end  of  each  record.   The  problem  was  due  to  internal  delays 
within  the  XDS  9300  which  caused  several  data  samples  to  be, 
omitted.   The  computer  laboratory  staff  revised  the  digiti- 
zation program  to  include  a  cross-check  between  each  sample 
and  lapsed  time. 

It  is  now  possible  to  sample  at  a  rate  of  at  least 
5000  SPS  and  check  each  block  for  missing  data.   Results 
obtained  from  PSD  analysis  indicate  the  problem  has  been 
corrected . 

3 .  PSD  Analysis  Procedures 

The  IBM  360/67  was  found  to  be  quite  capable  of  pro- 
cessing large  quantities  of  time  series  data.   The  operational 
routine  in  PSD  analysis  has  been  improved.  Methods  for 

affecting  faster  "turn  around"  have  been  developed. 

Several  methods  have  been  developed  for  statistically 
checking  the  digital  values  on  tapes  and  plotting  the  values 
by  the  Calcomp  plotter. 

4 .  PSD  Analysis  of  Turbulence  Signals 

A  definite  correlation  was  noted  between  the  temporal 
PSD  plots  of  both  temperature  and  velocity  and  the  original 
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analog  signal  from  which  the  PSD  results  were  computed.   This 
tended  to  further  support  the  conclusion  that  the  PSD  pro- 
grams were  working  properly. 

B.   RECOMMENDATIONS  FOR  FUTURE  WORK 

Due  to  the  fact  that  this  series  of  programs  has  proved  to 
be  an  extremely  powerful  tool  for  signal  analysis,  its  future 
use  should  be  vigourously  persued.   The  programs '  use  in  tur- 
bulence anlysis  has  been  well  established;  however,  its  ap- 
plication in  any  field  which  employs  PSD  techniques  should  be 
followed . 

The  time-varying  PSD  analysis  of  turblence  signal  should 
be  persued.   This  technique  showed  interesting  possibilities 
for  future  work.   Digital  tapes  of  four  turbulence  signals 
were  developed  and  can  be  easily  assessed  for  future  analysis 
of  the  signals.   Though  a  correlation  between  temperature  and 
velocity  was  noted,  new  digital  tapes  will  have  to  be  generated 
with  the  dual  channel  digitization  procedure  in  order  to  get 
cross-spectral  values  for  these  signals. 
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APPENDIX  C 


Equipment  Specifications 

1.   Gaussian  Noise  Generator 

Manufactured  by  Model  No.  Frequency  Ranges  Output  Spectrum 
Elgenco,  Inc.     603A       5Hz-20KHz        ±ldb 


(10Hz-500KHz) 


Output  Level     Spectral  Density  mv/  Hz 
20KHz  Setting     5  mv/  Hz 
3  V  81V    output  and 

20KH£mRange 


rms 


Filter 

Manufactured  by   Model  No.   Frequency  Range 
Khron-Hite,  Corp   33^0       .001  to  99-9  KHz 

Frequency  Accuracy   Pass  Band  Gain   Attenuation  Slope 
±2%  odb  or  20db      ' -48db/octave 
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